Gradient of image-space wave-equation tomography by the adjoint-state method

introduction

Wave-equation tomography aims to solve for earth models that explain observed seismograms under some norm. There are two main categories, depending on the domain in which the objective function is computed. In one category, known as waveform inversion (Woodward, 1992; Lines and Treitel, 1984; Tarantola, 1987), the objective function is defined in the data space, where the modeled data are compared with the recorded seismograms. In the other category, here called image-space wave-equation tomography (ISWET), the objective function is minimized in the image space. Two variants of ISWET are wave-equation migration velocity analysis (WEMVA) (Sava and Biondi, 2004a,b) and differential semblance velocity analysis (DSVA) (Shen, 2004; Shen and Symes, 2008).

Like waveform inversion, ISWET is a computationally demanding process. This computational cost is commonly decreased by using generalized sources (Tang et al., 2008; Shen and Symes, 2008). Because wavefield propagation is a linear process, generalized sources are computed by linearly combining source wavefields and receiver wavefields, using phase-encoding techniques (Romero et al., 2000; Whitmore, 1995). In particular, Guerra et al. (2009) used phase-encoding modeling in the image space to synthesize generalized source functions that drastically decrease the cost of DSVA.

WEMVA and DVSA seek the optimal slowness by driving an image perturbation to a minimum. However, they differ in the way the image perturbation is computed and, consequently, in the numerical optimization scheme applied. As Biondi (2008) points out, WEMVA is not easily automated. The image perturbation is computed by the linearized residual prestack-depth migration (Sava and Biondi, 2004a), which uses a manually picked residual-moveout parameter. Because the perturbed image computed with the linearized residual prestack-depth migration is consistent with the application of the forward wave-equation tomographic operator, the forward and adjoint ISWET operators can be used in conjugate-gradient methods to invert for the slowness perturbation. In DSVA, the perturbed image is computed by applying the fully automated differential-semblance operator (DSO) to the subsurface-offset gathers (ODCIGs) or angle gathers (ADCIGs). When applied to ODCIGs, DSVA minimizes the energy not focused at zero-offset. When applied to ADCIGs, DSVA minimizes energy departing from flatness of the reflectors. Although DSVA easily automates ISWET, it produces perturbed images that do not present the depth phase-shift introduced by the forward one-way ISWET operator, neither the DVSA amplitude behavior can be modeled by this operator. Therefore, the objective function computed with DSO is tipically minimized by quasi-Newton algorithms, which require computation of the gradient of the objective function.

The gradient of the objective function can be computed with the Frech $\acute {\rm e}$t derivatives. However, even for 2D applications of ISWET this computation can be very expensive. An efficient way to compute the gradient without using Frech $\acute {\rm e}$t derivatives is the adjoint-state method (Plessix, 2006; Chavent and Jacewitz, 1995). Plessix (2006) describes two methodologies for computing the gradient of the objective function using the adjoint-state method. One methodology uses the perturbation theory, which states that, at first order, a perturbation of the model parameters causes a perturbation of the objective function. The other uses the augmented Lagrangian. The augmented Lagrangian is formed by the objective function and the scalar product of the adjoint-state variables with general solutions of the forward modeling equations. The adjoint-state variables are, in turn, solutions of the adjoint-state equations. The adjoint-state equations are defined by equating to zero the derivatives of the augmented Lagrangian with respect to the state variables. For the linear case, the adjoint of the modeling operator applied to the adjoint-state variables gives the gradient of the objective function.

Although previous studies have computed the gradient of the ISWET objective function by the adjoint-state method (for example, Shen et al. (2003)), they have not provided a detailed derivation. Here, I show a detailed derivation the gradient of the ISWET objective function using the augmented Lagrangian methodology. The derivation is valid whether ISWET uses areal-shot migration or shot-profile migration. A complete description of the forward and adjoint of ISWET operators is given in Tang et al. (2008). A numerical example of slowness optimization on the Marmousi model illustrates the use of the gradient by a quasi-Newton algorithm.

Because I use image-space phase-encoded gathers in the numerical example, for completeness, I first briefly describe how to compute these phase-encoded gathers. A detailed treatment of the image-space phase-encoded gathers can be found in Biondi (2007,2006) and Guerra and Biondi (2008a). Then I derive the ISWET gradient using the adjoint-state method, and show the numerical example .

Gradient of image-space wave-equation tomography by the adjoint-state method