Wave-equation tomography in the generalized source domain

Yaxun Tang, Claudio Cardoso and Biondo Biondi September 11, 2008

Abstract

In this paper, we generalize wave-equation migration velocity analysis to composite shot domain. Any linear composite shot gather can be used for WEMVA. One particular case is the plane-wave shot gathers. I will demonstate some theoretical advantages of doing WEMVA in composite shot domain,

INTRODUCTION

WAVE-EQUATION TOMOGRAPHY

Wave-equation tomography is a non-linear inverse problem, it is based on an optimization that tries to find an optimal background slowness that minimizes the residual field, 1I , defined in the image space. The residual field is derived from the background image, I, which is computed with a background slowness (or the current estimate of the slowness). The residual field measures the correctness of the background slowness, its minimum (under some norm, e.g. ℓ_{2}) is achieved when a correct background slowness has been used for migration. There are many choices of the residual field, such as residual moveout in the Angle-Domain Common-Image Gathers (ADCIGs), differential semblence in the ADCIGs, reflection angle stacking power (in this case we have to maximize the residual field, or minimize the minus stacking power), and etc.. For simplicity, we define a linear operator K that measures the image perturbation as follows:

1I D (1 − K)I,  (1)

where 1 is the identity operator. In the wave-equation migration velocity analysis (WEMVA) method, K is the differential residual migration operator (Sava and Biondi???), it flattens the ADCIGs and also focuses the migrated image in the CMP direction; In the Differential Semblance Optimization (DSO) method,

K D 1 − R,  (2)

where R is the DSO operator either in the subsurface offset domain or in the angle domain. Subsurface offset domain DSO focuses the energy at zero offset while angle domain DSO flattens the ADCIGs.

In general, if we choose ℓ_{2} norm, the objective function to minimize can be written as follows:

J D \frac{1}{2}||1I ||_{2} D \frac{1}{2}||(1 − K)I||_{2},  (3)

where || · ||_{2} stands for the ℓ_{2} norm. Gradient based optimization techniques such as Quasi-Newton method and Conjugate Gradient method can be used to minimize the objective function J, which is non-linear because the image I is a non-linear function of the slowness s. To use the gradient based methods, we have to provide the gradient of J with respect to the slowness s. The gradient reads

∇J D ℜ( ( \frac{∂I}{∂s}) ′(1 − K)′(1 − K)I) ,  (4)

where′means the adjoint, \frac{∂I}{∂s} is a linear operator which defines the mapping from the slowness perturbation 1s to the image perturbation 1I . This is easy to see by expanding the image I around a background slowness ŝ as follows:

I D ÎC \frac{∂I}{∂s}(s −ŝ) C ··· ,  (5)

where Î is the background image computed with the background slowness ŝ. Keeping only the zero and first order terms, we get the linearized operator \frac{∂I}{∂s} defined as follows:

1I D \frac{∂I}{∂s}1s D T1s ,  (6)

where 1I D I −Î and 1s D s −ŝ. T D \frac{∂I}{∂s} is the wave-equation tomographic operator. The tomographic operator can either be defined in the source and receiver domain (Sava,????) or in the shot-profile domain (Shen,???). In next section we follow a similar approach disscussed by Shen,?? and review the forward and adjoint tomographic operator in the shot profile domain, in the subsequent sections, we generalize the expression of the tomographic operator to generalized source domains.

THE TOMOGRAPHIC OPERATOR IN THE SHOT-PROFILE DOMAIN

For the conventional shot-profile migration, both source and receiver wavefields are downward continued with the following one-way wave equations (Claerbout,???):

\{ \begin{matrix}( \frac{∂}{∂z} − \sqrt{ω^{2}s^{2}(x) − |k|^{2}} ) D(x,x_{s},ω) D 0 \\ D(x,y,z D 0,x_{s},ω) D f_{s}(ω)δ(x − x_{s}) \end{matrix} ,  (7)

and

\{ \begin{matrix}( \frac{∂}{∂z} C \sqrt{ω^{2}s^{2}(x) − |k|^{2}} ) U(x,x_{s},ω) D 0 \\ U(x,y,z D 0,x_{s},ω) D Q(x,y,z D 0,x_{s},ω) \end{matrix} ,  (8)

where D(x,x_{s},ω) is the source wavefield for a single frequency ω at image point x D (x,y,z) with the source located at x_{s} D (x_{s},y_{s},0), U(x,x_{s},ω) is the receiver wavefield for a single frequency ω at image point x for the source located at x_{s}. s(x) is the slowness at x, k D (k_{x},k_{y}) is the spatial wavenumber vector, f_{s}(ω) is the freqnency dependent source signature and f_{s}(ω)δ(x − x_{s}) defines the point source function at x_{s}, which serves as the boundary condition of equation 7. Q(x,y,z D 0,x_{s},ω) is the shot gather for the shot located at x_{s}, which serves as the boundary condition of equation 8. To produce the image, the following cross-correlation imaging condition is used:

I(x,h) D ∑_{x_{s}} ∑_{ω} \overline{D(x − h,x_{s},ω)} U(x C h,x_{s},ω),  (9)

where h D (h_{x},h_{y},h_{z}) is the subsurface half offset.

The perturbed image at some image point can be derived by a simple application of the chain rule to equation 9:

1I(x,h) D ∑_{x_{s}} ∑_{ω} ( \overline{1D(x − h,x_{s},ω)} Û(x C h,x_{s},ω)C
\overline{D̂(x − h,x_{s},ω)} 1U(x C h,x_{s},ω)) ,  (10)

where D̂(x − h,x_{s},ω) and Û(x C h,x_{s},ω) are the background source and receiver wavefields computed with the background slowness ŝ(x); 1D(x − h,x_{s},ω) and 1U(x C h,x_{s},ω) are the perturbed source wavefield and perturbed receiver wavefield, which are the results of the slowness perturbation 1s(x). The perturbed source and receiver wavefields satisfy the following linearized (with respect to slowness) one-way wave equations, respectively(see Appendix A for derivations)

\{ \begin{matrix}( \frac{∂}{∂z} − \sqrt{ω^{2}s^{2}(x) − |k|^{2}} ) 1D(x,x_{s},ω) D _{\frac{Ciω1s(x)}{\sqrt{1− \frac{|k|^{2}}{ω^{2}ŝ(x)}}}} D̂(x,x_{s},ω) \\ 1D(x,y,z D 0,x_{s},ω) D 0 \end{matrix} ,  (11)

and

\{ \begin{matrix}( \frac{∂}{∂z} C \sqrt{ω^{2}s^{2}(x) − |k|^{2}} ) 1U(x,x_{s},ω) D _{\frac{−iω1s(x)}{\sqrt{1− \frac{|k|^{2}}{ω^{2}ŝ(x)}}}} Û(x,x_{s},ω) \\ 1U(x,y,z D 0,x_{s},ω) D 0 \end{matrix} .  (12)

Equations 11 and 12 are linear with respect to slowness perturbation, thus they define linear mappings from the slowness perturbation 1s(x) to the perturbed source wavefield 1D(x,x_{s},ω) and the perturbed receiver wavefield 1U(x,x_{s},ω), respectively. Therefore, recursively solving equations 11 and 12 give us the perturbed source and receiver wavefields. The perturbed source and receiver wavefields are then used in equation 10 to generate the perturbed image I(x,h), where the background source and receiver wavefields are precomputed by solving the recursive equations 7 and 8 with a background slowness ŝ(x). Appendix B gives a more detailed matrix representation of how to evaluate the forward tomographic operator T.

To compute the adjoint tomographic operator T′, we first apply the adjoint of the imaging condition in equation 10 to get the perturbed source and receiver wavefields 1D(x,x_{s},ω) and 1U(x,x_{s},ω) as follows

1D(x,x_{s},ω) D ∑_{h} 1I(x,h)Û(x C h,x_{s},ω)  (13)
1U(x,x_{s},ω) D ∑_{h} 1I(x,h)D̂(x − h,x_{s},ω)  (14)

Then we solve the adjoint equations of 11 and 12 to get the slowness perturbation 1s(x). Again, in order to solve the adjoint equations of 11 and 12, the background source wavefield D̂(x,x_{s},ω) and the background receiver wavefield Û(x,x_{s},ω) have to be computed in advance. Appendix B gives a more detailed matrix representation of how to evaluate the adjoint tomographic operator: T′.

TOMOGRAPHY WITH THE ENCODED WAVEFIELDS

However, the cost for computing the gradient of the objective function J in the shot-profile domain is twice the cost of a shot-profile migration, to minimize J many gradient and function evaluations have to be computed, thus wave-equation tomography in the conventional shotprofile domain seems infeasible with nowadays computer resources. To save the cost and to make this powerful method more practical, we extend the wave-eqaution tomography theory to generalized source domains, or phase encoded source domain. We discuss two different strategies generate the generalized shot gathers, i.e.,

Data-space encoded wavefields

We encode the source wavefield and the corresponding receiver wavefield

D̃(x,p_{s},ω) D ∑_{x_{s}} D(x,x_{s},ω)α(x_{s},p_{s},ω)  (15)
Ũ(x,p_{s},ω) D ∑_{x_{s}} U(x,x_{s},ω)α(x_{s},p_{s},ω)  (16)

where α(x_{s},p_{s},ω) is some phase encoding function. In the case of plane-wave phase encoding, α(x_{s},p_{s},ω) is defined as

α(x_{s},p_{s},ω) D e^{iωp_{s}x_{s}} ,  (17)

where p_{s} is the ray parameter for the source plane waves on the surface. Because of the linearity of the wave equation, the en

Imaging in the generallized source domain

I(x,h) D ∑_{p_{s}} ∑_{ω} D̃(x − h,p_{s},ω)Ũ(x C h,p_{s},ω)  (18)

Image-space encoded wavefields

NUMERICAL EXAMPLES

This appendix derives the linearized one-way wave equation

This appendix demonstrates a matrix representation of the forward tomographic operator T

D_{zC1z} D E_{z}(s_{z})1D _{z},  (19)

where

E_{z}(s_{z}) D e^{ik_{z}(s_{z})1z} D e^{i\sqrt{ω^{2}s_{z}^{2}−|k|^{2}}}  (20)

The perturbed source wavefield at some depth level can be derived from the background wavefield by a simple application of the chain rule to equation 19:

1D _{zC1z} D E_{z}(ŝ_{z})1D _{z}C 1E _{z}(ŝ_{z})D_{z},  (21)

where 1E _{z} represents the perturbed extrapolator, it can be obtained by a formal linearization with respect to slowness of the extrapolator defined in equation 20

E_{z}(s_{z}) D e^{ik_{z}(s_{z})1z} ≈ e^{i1zk̂_{z}} C e^{i1zk̂_{z}} ( i1z \frac{dk_{z}}{ds_{z}}|_{s_{z}Dŝ_{z}}) 1s _{z}

D E_{z}(ŝ_{z}) C E_{z}(ŝ_{z})( i1z \frac{dk_{z}}{ds_{z}}|_{s_{z}Dŝ_{z}}) 1s _{z}  (22)

where k̂_{z} D k_{z}(ŝ_{z}). So the perturbed extrapolator reads

1E _{z}(ŝ_{z}) D E_{z}(ŝ_{z})( i1z \frac{dk_{z}}{ds_{z}}|_{s_{z}Dŝ_{z}}) 1s _{z}  (23)

Substitute equation 23 into 21, we get

1D _{zC1z} D E_{z}(ŝ_{z})1D _{z}C E_{z}(ŝ_{z})( i1z \frac{dk_{z}}{ds_{z}}|_{s_{z}Dŝ_{z}}) D_{z}1s _{z},  (24)

Let us define a scattering operator G_{z} that interacts with the background wavefield as follows:

G_{z}(D_{z},ŝ_{z}) D ( i1z \frac{dk_{z}}{ds_{z}}|_{s_{z}Dŝ_{z}}) D_{z} D \frac{iω1z}{\sqrt{1 − \frac{|k|^{2}}{ω^{2}ŝ_{z}^{2}}}} D_{z}  (25)

Then the perturbed source wavefield for depth level z C 1z can be rewritten as follows:

1D _{zC1z} D E_{z}(ŝ_{z})1D _{z}C E_{z}(ŝ_{z})G_{z}(D̂_{z},ŝ_{z})1s _{z}.  (26)

We can further write out the recursive equation 26 for all depth levels in the following matrix form:

( 1D _{0} ; 1D _{1} ; 1D _{2} ; ⋮ ; 1D _{n} ) D ( \begin{matrix}0 & 0 & 0 & ··· & 0 & 0 \\ E_{0} & 0 & 0 & ··· & 0 & 0 \\ 0 & E_{1} & 0 & ··· & 0 & 0 \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ \\ 0 & 0 & 0 & ··· & E_{n−1} & 0 \end{matrix} ) ( 1D _{0} ; 1D _{1} ; 1D _{2} ; ⋮ ; 1D _{n} ) C

( \begin{matrix}0 & 0 & 0 & ··· & 0 & 0 \\ E_{0} & 0 & 0 & ··· & 0 & 0 \\ 0 & E_{1} & 0 & ··· & 0 & 0 \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & 0 \\ 0 & 0 & 0 & ··· & E_{n−1} & 0 \end{matrix} ) ( \begin{matrix}G_{0} & 0 & 0 & ··· & 0 \\ 0 & G_{1} & 0 & ··· & 0 \\ 0 & 0 & G_{2} & ··· & 0 \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ \\ 0 & 0 & 0 & ··· & G_{n} \end{matrix} ) ( 1s _{0} ; 1s _{1} ; 1s _{2} ; ⋮ ; 1s _{n} )

Or in a more compact notation

1D D E(ŝ)1D C E(ŝ)G(D̂,ŝ)1s .  (27)

The solution of equation 27 can be formally written as follows:

1D D (1 − E(ŝ)) ^{−1} E(ŝ)G(D̂,ŝ)1s .  (28)

Similarly, the receiver wavefield satisfy the following recursive relation:

U_{zC1z} D \overline{E_{z}(ŝ_{z})}1U _{z},  (29)

where the over line stands for the complex conjugate, thus

\overline{E_{z}(s_{z})} D e^{−ik_{z}(s_{z})1z} D e^{−i\sqrt{ω^{2}s_{z}^{2}−|k|^{2}} 1z} .  (30)

Apply the chain rule to euqation 29, we obtain the recursive equation for the scattered receiver wavefield:

1U _{zC1z} D \overline{E_{z}(ŝ_{z})}1U _{z}− \overline{E_{z}(ŝ_{z})}G_{z}(Û_{z},ŝ_{z})1s _{z},  (31)

where G_{z}(Û_{z},ŝ_{z}) is the scattering operator interacts with the background receiver wavefield as follows:

G_{z}(Û_{z},ŝ_{z}) D ( i1z \frac{dk_{z}}{ds_{z}}|_{s_{z}Dŝ_{z}}) Û_{z} D \frac{iω1z}{\sqrt{1 − \frac{|k|^{2}}{ω^{2}ŝ_{z}^{2}}}} Û_{z}.  (32)

We can further write out the recursive equation 32 for all depth levels in the following matrix form:

( 1U _{0} ; 1U _{1} ; 1U _{2} ; ⋮ ; 1U _{n} ) D ( \begin{matrix}0 & 0 & 0 & ··· & 0 & 0 \\ \overline{E_{0}} & 0 & 0 & ··· & 0 & 0 \\ 0 & \overline{E_{1}} & 0 & ··· & 0 & 0 \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ \\ 0 & 0 & 0 & ··· & \overline{E_{n−1}} & 0 \end{matrix} ) ( 1U _{0} ; 1U _{1} ; 1U _{2} ; ⋮ ; 1U _{n} ) −

( \begin{matrix}0 & 0 & 0 & ··· & 0 & 0 \\ \overline{E_{0}} & 0 & 0 & ··· & 0 & 0 \\ 0 & \overline{E_{1}} & 0 & ··· & 0 & 0 \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & 0 \\ 0 & 0 & 0 & ··· & \overline{E_{n−1}} & 0 \end{matrix} ) ( \begin{matrix}G_{0} & 0 & 0 & ··· & 0 \\ 0 & G_{1} & 0 & ··· & 0 \\ 0 & 0 & G_{2} & ··· & 0 \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ \\ 0 & 0 & 0 & ··· & G_{n} \end{matrix} ) ( 1s _{0} ; 1s _{1} ; 1s _{2} ; ⋮ ; 1s _{n} )

Or in a more compact notation

1U D \overline{E(ŝ)}1U − \overline{E(ŝ)}G(Û,ŝ)1s .  (33)

The solution of equation 33 can be formally written as follows:

1U D (1 − \overline{E(ŝ)}) ^{−1} \overline{E(ŝ)}G(Û,ŝ)1s .  (34)

Therefore, the perturbed image can be obtained as follows:

( 1I _{0} ; 1I _{1} ; 1I _{2} ; ⋮ ; 1I _{n} ) D ( \begin{matrix}Û_{0} & 0 & 0 & ··· & 0 \\ 0 & Û_{1} & 0 & ··· & 0 \\ 0 & 0 & Û_{2} & ··· & 0 \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ \\ 0 & 0 & 0 & ··· & Û_{n} \end{matrix} ) ( \overline{1D _{0}} ; \overline{1D _{1}} ; \overline{1D _{2}} ; ⋮ ; \overline{1D _{n}} ) C ( \begin{matrix}\overline{D̂}_{0} & 0 & 0 & ··· & 0 \\ 0 & \overline{D̂_{1}} & 0 & ··· & 0 \\ 0 & 0 & \overline{D̂_{2}} & ··· & 0 \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ \\ 0 & 0 & 0 & ··· & \overline{D̂_{n}} \end{matrix} ) ( 1U _{0} ; 1U _{1} ; 1U _{2} ; ⋮ ; 1U _{n} )

Or in a more compact notation

1I D diag\binom{}{U}\overline{1D} C \overline{diag\binom{}{D}}1U .  (35)

After substituting equations 28 and 34 into equation 35

1I D (diag\binom{}{U}\overline{(1 − E(ŝ)) ^{−1} E(ŝ)G(D̂,ŝ)} C
\overline{diag\binom{}{D}}(1 − \overline{E(ŝ)}) ^{−1} \overline{E(ŝ)}G(Û,ŝ)) 1s ,  (36)

we finally get the following matrix representation of the forward tomographic operator T:

T D diag\binom{}{U}\overline{(1 − E(ŝ)) ^{−1} E(ŝ)G(D̂,ŝ)} C
\overline{diag\binom{}{D}}(1 − \overline{E(ŝ)}) ^{−1} \overline{E(ŝ)}G(Û,ŝ)  (37)

This appendix demonstrates a matrix representation of the adjoint tomographic operator T′ It is easy to get the adjoint equation from equation 36,