Wave-equation tomography in the generalized source domain

Yaxun Tang, Claudio Guerra and Biondo Biondi

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 (?), 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 (?) or in the shot-profile domain (?). In next section we follow a similar approach disscussed by (?) 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} − i\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 i\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} − i\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 i\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

As introduced by (?), the prestack exploding reflector modeling aims to generate a smaller dataset to speed up iterations of migration velocity analysis. The method synthesizes areal source and receiver data at the surface by recursively upward propagating source and receiver wavefields, having prestack images obtained by wave-equation migration with an inaccurate velocity, ŝ(x), as initial conditions.

The reduction of the dataset size to be migrated is achieved by modeling several wavefields at the same time and collecting them into super-areal data. To avoid crosstalk, due to correlation of unrelated wavefields during migration, phase encoding is applied in the modeling step. (?) showed that random phases can be used to encode the modeling experiments such that the crosstalk is dispersed throughout the image. By performing more modeling realizations, the signal to noise ratio is improved.

For the areal shot migration of the synthesized data, both synthesized source and receiver wavefields are downward continued with the following one-way wave equations:

\{ \begin{matrix}( \frac{∂}{∂z} − i\sqrt{ω^{2}ŝ^{2}(x) − |k|^{2}} ) D̂(x,p_{m},ω) D 0 \\ D̂(x,y,z D 0,p_{m},ω) D D̆(x,y,z D 0,p_{m},ω) \end{matrix} ,  (19)

and

\{ \begin{matrix}( \frac{∂}{∂z} C i\sqrt{ω^{2}ŝ^{2}(x) − |k|^{2}} ) Û(x,p_{m},ω) D 0 \\ Û(x,y,z D 0,p_{m},ω) D Ŭ(x,y,z D 0,p_{m},ω) \end{matrix} ,  (20)

where D̆ and Ŭ defines, respectively, modeling realizations of the synthesized source and receiver wavefields at the surface. They are used as boundary conditions for the downward continuation. D̂ and Û are the downward propagated background source and receiver wavefields, respectively. The vector p_{m} represents the modeling realizations.

The background image is computed by applying the cross-correlation imaging condition according to:

Î(x,h) D ∑_{p_{m}} ∑_{ω} \overline{D̂(x − h,p_{m},ω)} Û(x C h,p_{m},ω).  (21)

Ideally, if the phase encoding scheme is effective, the image Î(x,h) shows no crosstalk.

Similarly, to the shot-profile case, a perturbed image, 1I, can be obtained by applying the chain rule to equation ??:

1I(x,h) D ∑_{p_{m}} ∑_{ω} ( \overline{1D(x − h,p_{m},ω)} Û(x C h,p_{m},ω)C
\overline{D̂(x − h,p_{m},ω)} 1U(x C h,p_{m},ω)) ,  (22)

where the wavefields indicated with a ˆ . are the background wavefields and the perturbed source and receiver wavefields are 1D and 1U, respectively. The perturbed wavefields satisfy the one-way wave equation:

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

and

\{ \begin{matrix}( \frac{∂}{∂z} − i\sqrt{ω^{2}ŝ^{2}(x) − |k|^{2}} ) 1U (x,p_{m},ω) D _{\frac{−iω1ŝ(x)}{\sqrt{1− \frac{|k|^{2}}{ω^{2}ŝ^{2}(x)}}}} Û(x,p_{m},ω) \\ 1D (x,y,z D 0,p_{m},ω) D 0 \end{matrix} ,  (24)

D(x,ω) D [D(x,x_{m},ω) C I_{D}(x_{m},h)] e^{−ik_{z}1z}  (25)
U(x,ω) D [U(x,x_{m},ω) C I_{U}(x_{m},h)] e^{ik_{z}1z}.  (26)

Here The different initial conditions, I_{D} and I_{U}, for the source and receiver wavefields respectively, are computed by a proper rotation of the original prestack image necessary for the case of using an inaccurate velocity and in the presence of dipping layers (?).

Reduction of the dataset size to be migrated is achieved by modeling several wavefields at the same time and collecting them into a super-areal data. The wavefields resulting from the modeling realizations are computed by

D̃(x,p_{m},ω) D ∑_{x_{m}} [D(x,z C 1z,ω) C α(x_{m},p_{m},ω)I_{d}(x − h,z C 1z)] e^{−ik_{z}1z} (27)
Ũ(x,p_{m},ω) D ∑_{x_{m}} [U(x,z C 1z,ω) C α(x_{m},p_{m},ω)I_{u}(x C h,z C 1z)] e^{ik_{z}1z}  (28)

where α(x_{m},p_{m},ω) D e^{iγ (x_{m},p_{m},ω)} is the phase encoding function, x_{m} denotes model space coordinates and p_{m} represents the realizations. γ is any function that minimizes the crosstalk. In the present case, we use a pseudo-random function following (?).

The drawback of the data reduction is crosstalk generation when applying the imaging condition during migration. To attenuate crosstalk, source and receiver wavefields are encoded during the modeling step. (?) showed that reasonable attenuation can be achieved by applying random phase encoding in (x,z,ω)-dimensions. Furthermore, by migrating more modeling realizations, the signal to noise ratio is improved.

NUMERICAL EXAMPLES

REFERENCES