Gradient of image-space wave-equation tomography in the generalized

source domain by the adjoint-state method

Claudio Guerra Stanford Exploration Project, Mitchell Bldg., Department of Geophysics, Stanford University, Stanford, CA 94305-2215 (September 24, 2009)

ABSTRACT

Optimization using gradient descent techniques requires the computation of the gradient of the objective function. The adjoint-state method is an efficient way of computing the gradient, since is not necessary to compute the expensive Frechét derivatives. Here, I derive the gradient of the image-space wave-equation tomography using the adjoint state method and show its application with a numerical example.

INTRODUCTION

Wave-equation tomography aims to solve for Earth models that explain the data under some norm. In general, two categories are distinguished depending on the domain the objective function is computed. In one category, known as waveform inversion (???), the objective function is defined in the data space where the modeled data are compared with the recorded seismograms. In the other category, generally called here as image-space wave-equation tomography (ISWET) from which waveequation migration velocity analysis (WEMVA) (??) and differential semblance velocity analysis (DSVA) (??) are two variants, the objective function to be minimized is defined in the image space.

Despite their differences, those two categories are linked by the concept of extended modeling (?). ? showed that ISWET can be regarded as a solution method for the partially linearized waveform inversion. While ISWET solves for the long-wavelength components of the slowness model, assuring global convergence, waveform inversion computes the short-wavelength components.

Similar to waveform inversion, ISWET is a computationally demanding process. To decrease its computational cost it is common practice to use generalized sources (??). Because wavefield propagation is a linear process, generalized sources can be computed by linearly combining source wavefields and receiver wavefields, using phase encoding techniques (??). In particular, ? used generalized source functions synthesized by phase encoding modeling experiments in the image space drastically decreasing the cost of DSVA. The example showed in this report also uses this kind of generalized source functions.

The aforementioned ISWET methods differ on the way the objective function is computed and, consequently, the numerical optimization scheme applied. As pointed out by ?, in WEMVA, an intermediate residual-moveout parameter picking is necessary to produce the perturbed image via linearized residual prestack depth migration (?). The picking procedure makes it difficult to automate the whole process. As the computed perturbed image is consistent with the wave-equation tomographic operator, however, conjugate gradient methods can be used to invert for the slowness perturbation. In DSVA, the perturbed image is computed by applying the fully automated differential-semblance operator (DSO) on the background image, both using the subsurface-offset gathers (ODCIG) or angle gathers (ADCIG). When applied to ODCIG, DSVA minimizes the energy not focused at zerooffset. When applied to ADCIG, DSVA minimizes energy departing from flatness of the reflectors. Notwithstanding the ability of automating ISWET, DSVA produces perturbed images not consistent with the wave-equation tomographic operator. To minimize the objective function computed with DSO, quasi-Newton algorithms can be used for which the gradient of the objective function needs to be computed.

The gradient of the objective function can be obtained by computing the Frechét derivatives. However, this computation even for 2D applications of ISWET can be very expensive. An efficient way of computing the gradient without the need of Frechét derivatives is the adjoint-state method (??). ? describes two methodologies for computing the gradient of the objective function using the adjoint-state method. The more intuitive is the one that uses the augmented Lagrangian. The augmented Lagrangian is formed by the objective function and the scalar product of the adjointstate variables and general solutions of the forward modeling equations. The adjoint-state variables are solution of a system of adjoint-state equations defined by equating to zero the derivative 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.

Here, I show how to compute the gradient of the ISWET objective function using the augmented Lagrangian methodology. The derivation of the gradient presented here is valid whether ISWET uses areal-shot migration or shot-profile migration. ? provide a complete derivation of the forward and adjoint wave-equation tomographic operators in the generalized source domain. The goal of this paper is to provide a more detailed derivation of the gradient of the ISWET objective function using the adjoint state method than the ones available, for example, ?.

The references on the image-space phase-encoded gathers are ?? and ?. For completeness, here I briefly describe how to compute these phase-encoded gathers. Then, I derive the ISWET gradient using the adjoint-state method and, finally show an example of DSVA using image-space phaseencoded gathers on the Marmousi model.

IMAGE-SPACE GENERALIZED SOURCES

The prestack exploding-reflector modeling, introduced by ? and ?, synthesizes areal data and the corresponding areal source function, having as initial condition, a wave-equation migrated prestack image obtained with an inaccurate slowness model. In general, one single reflection event from one single ODCIG is modeled by the recursive upward continuation with the following one-way wave equations:

\{ \begin{matrix}( \frac{∂}{∂z} − i\sqrt{ω^{2}ŝ^{2}(x) − |k|^{2}} ) d(x,ω;x_{m},y_{m}) D r_{D}(x,h;x_{m},y_{m}) \\ d(x,y,z D z_{max},ω;x_{m},y_{m}) D 0 \end{matrix} ,  (1)

and

\{ \begin{matrix}( \frac{∂}{∂z} C i\sqrt{ω^{2}ŝ^{2}(x) − |k|^{2}} ) u(x,ω;x_{m},y_{m}) D r_{U}(x,h;x_{m},y_{m}) \\ u(x,y,z D z_{max},ω;x_{m},y_{m}) D 0 \end{matrix} ,  (2)

where r_{D}(x,h;x_{m},y_{m}) and r_{U}(x,h;x_{m},y_{m}) are the isolated SODCIGs at the horizontal location (x_{m},y_{m}) for a single reflector, and are suitable for the initial conditions for the source and receiver wavefields, respectively. They are obtained by rotating the original unfocused SODCIGs according to the apparent geological dip of the reflector. This rotation maintains the velocity information needed for migration velocity analysis, especially for dipping reflectors (?). d(x,y,z D z_{c},ω;x_{m},y_{m}) is the areal source data and u(x,y,z D z_{c},ω;x_{m},y_{m}) is the areal receiver data for a single reflector and a single SODCIG located at (x_{m},y_{m}). z D z_{c} denotes that the wavefields can be collected at any depth level, z_{c}. This characteristic is important to accelerate ISWET especially if z_{c} is such that it separates regions of sufficiently accurate slowness above and inaccurate slowness below. As a consequence, the synthesized gathers are naturally ‘datumized’ and the wavefield propagations during ISWET can be restricted to the region where the slowness model must be updated. This feature allows ISWET be applied in a target-oriented manner.

However, if one considers to model one single reflection from one single SODCIG at a time, the modeled dataset can be orders of magnitude bigger than the original data set. As discussed by ???, by modeling several reflectors and several SODCIGs simultaneously using random phase encoding a much smaller dataset is obtained. The randomly encoded areal source and areal receiver wavefields can be computed as follows:

\{ \begin{matrix}( \frac{∂}{∂z} − i\sqrt{ω^{2}ŝ^{2}(x) − |k|^{2}} ) d̃(x,p_{m},ω) D r̃_{D}(x,h,p_{m},ω) \\ d̃(x,y,z D z_{max},p_{m},ω) D 0 \end{matrix} ,  (3)

and

\{ \begin{matrix}( \frac{∂}{∂z} C i\sqrt{ω^{2}ŝ^{2}(x) − |k|^{2}} ) ũ(x,p_{m},ω) D r̃_{U}(x,h,p_{m},ω) \\ ũ(x,y,z D z_{max},p_{m},ω) D 0 \end{matrix} ,  (4)

where r̃_{D}(x,h,p_{m},ω) and r̃_{U}(x,h,p_{m},ω) are the encoded SODCIGs after rotations. They are defined as follows:

r̃_{D}(x,h,p_{m},ω) D ∑_{x_{m}} ∑_{y_{m}} r_{D}(x,h,x_{m},y_{m})β(x,x_{m},y_{m},p_{m},ω),  (5)
r̃_{U}(x,h,p_{m},ω) D ∑_{x_{m}} ∑_{y_{m}} r_{U}(x,h,x_{m},y_{m})β(x,x_{m},y_{m},p_{m},ω),  (6)

where β(x,x_{m},y_{m},p_{m},ω) D e^{iγ (x,x_{m},y_{m},p_{m},ω)} is chosen to be the random phase-encoding function, with γ (x,x_{m},y_{m},p_{m},ω) being a uniformly distributed random sequence in x, x_{m}, y_{m} and ω; the variable p_{m} is the index of different realizations of the random sequence. Recursively solving Equations ?? and ?? gives us the encoded areal source data d̃(x,y,z D z_{c},p_{m},ω) and areal receiver data ũ(x,y,z D z_{c},p_{m},ω), which can be collected at the depth level, z_{c}.

GRADIENT COMPUTATION OF ISWET BY THE ADJOINT STATE METHOD

Image-space wave-equation tomography aims to solve for the slowness model, s D s(x), that minimizes the linearized objective function

J (s) D \frac{1}{2} ||1r(s)| |^{2} D \frac{1}{2} ||r(s) − Mr(s)| |^{2},  (7)

where 1r D 1r(x,h) is the image perturbation which measures the goodness of the slowness model. Notice that the objective function depends on the slowness model through the image perturbation. 1r is computed by applying a differential residual-focusing operator M (?), using either differential residual prestack-migration (??) or differential-semblance optimization (DSO) (?) operators to the image r D r(x,h). Herein, operators are represented by bold capital letters.

If the differential residual-focusing operator M is independent on the slowness, the gradient of this objective function, evaluated at the current slowness, ŝ D ŝ(x), is

∇J (s) D ( \frac{∂r}{∂s}) ′|_{sDŝ} (I − M ′)1r̂ .  (8)

where ‘′’ denotes the adjoint, I is the identity operator and 1r̂ D 1r̂(x,h) is the perturbed image obtained with the current slowness model. The linear operator \frac{∂r̂}{∂s} defines the mapping, \frac{∂r̂}{∂s}1s D 1r, between the slowness perturbation 1s and the image perturbation and it is called image-space waveequation tomographic operator.

As the image-space wave-equation tomographic operator is composed of different operators, it is difficult from equation ?? to envision which operations are performed to compute the gradient. Therefore, for a clear explanation of the operators involved, I use the adjoint-state method to derive the gradient of the objective function (equation ??).

In migration with generalized sources (areal shot and shot-profile migration, e.g.), the source and receiver wavefields are propagated independently and the image, r_{z} D r_{z}(x,h), at a depth level z, is computed by the crosscorrelation

r_{z}(x,h) D ∑_{ω} d_{z}^{∗}(x − h,ω)u_{z}(x C h,ω),  (9)

where d_{z}(x,ω) is the source wavefield for a single frequency ω at horizontal coordinates x D (x,y) ; u_{z}(x,ω) is the receiver wavefield and h D (h_{x},h_{y}) is the subsurface half-offset, and ‘^{∗}’ stands for the complex-conjugate. An additional summation over shots is required when migrating more than one shot. Hereafter, letters d and u stand for source and receiver wavefields, respectively, irrespective to the migration scheme.

In a more compact notation, not explicitly writing the dependencies on x and h, equation ?? can be written as:

r_{z} D SD′_{z}(ω)u_{z}(ω) D SU_{z}(ω)d_{z}^{∗}(ω),  (10)

where D and U are convolutional operators composed of (h_{x},h_{y})-shifted versions of d_{z}(x,ω) and u_{z}(x,ω), respectively. Operator S corresponds to the summation over frequency.

For subsequent depth levels, d(x,ω) is computed by means of the recursive downward propagation

\{ \begin{matrix}d_{zC1}(ω) D T_{z}^{↓}(ω,s)d_{z}(ω) \\ d_{1}(ω) D q(ω), \end{matrix}  (11)

where T_{z}^{↓} is the downward continuation operator, which is function of the slowness, s, and q(ω) is the source wavefield used as boundary condition. In the case of conventional shot-profile migration q(ω) D f_{s}(ω)δ(x − x_{s} ) is the source signature located at x_{s} D (x_{s},y_{s},0). If using the generalized sources in the image space, q(ω) represents the image-space phase-encoded source wavefield of equation ??.

The downward continuation of the receiver wavefield is performed by

\{ \begin{matrix}u_{zC1}(ω) D T_{z}^{↓}(ω,s)u_{z}(ω) \\ u_{1}(ω) D w(ω), \end{matrix}  (12)

where w(ω) is the recorded data at the surface for shot-profile migration. If using generalized sources in the image space, w(ω) is the phase-encoded areal receiver wavefield of equation??. In equations ?? and ??, I omitted the dependencies of the wavefield with respect to x. The subscript 1 in equations ?? and ?? represents the surface for the shot-profile migration and the ‘collection’ depth level, z_{c}, for the image-space phase-encoded wavefields.

In the image-space wave-equation tomography problem, the perturbed source and receiver wavefields and image perturbations are used to compute the slowness perturbation to update the current slowness model. From the perturbation theory, we have that d D d̂C 1d, u D ûC 1u and, consequently, r D r̂ C1r are physical realizations with s D ŝC1s, where the hat refers to fields obtained with the background slowness. To the first order (Born approximation), these perturbed fields are given by

1d_{zC1}(ω) D T_{z}^{↓}(ω,ŝ)1d_{z}(ω) C˜D̃_{z}(ω)1s_{z}  (13)
1u_{zC1}(ω) D T_{z}^{↓}(ω,ŝ)1u_{z}(ω) C˜Ũ_{z}(ω)1s_{z}  (14)

The diagonal operators ˜ D̃_{z} and ˜ Ũ_{z} have in the diagonal entries the scattered source and receiver wavefields, respectively, are given by the action of the scattering operator 1T_{z}^{↓} on the background wavefields

˜D̃_{z}(ω) D 1T_{z}^{↓}(ω,ŝ)d̂_{z}(ω) D i \frac{ω^{2}ŝ}{\sqrt{ω^{2}s^{2}− |k|^{2}}} dz d̂_{z}(ω),  (15)

and

˜Ũ_{z}(ω) D 1T_{z}^{↓}(ω,ŝ)û_{z}(ω) D −i \frac{ω^{2}ŝ}{\sqrt{ω^{2}s^{2}− |k|^{2}}} dz û_{z}(ω).  (16)