RESEARCH PROPOSAL Wave equation velocity inversion using exploding reflectors modeling

Claudio Guerra 10:00am, Monday June 16, 2008

Ph.D Proposal  2  Guerra

ABSTRACT

Problem

To obtain a seismic image of the Earth’s subsurface, estimates of the seismic velocity are required at different steps of seismic processing. The accuracy of the results depends on how precisely the velocity estimates represent the subsurface, based on the physical assumptions of the wave propagation.

Velocity model estimation involves some transformation of the data and interpretation of velocity parameters. For instance, the normal-moveout (NMO) velocity is obtained by fitting hyperbolas to data sorted to common-midpoint (CMP). Interval velocity for depth migration is determined iteratively by: 1) transforming data into common-image gathers (CIG) by means of prestack-depth migration with a presumed velocity model; 2) measuring the misfocusing on CIG; and 3) updating the velocity model based on the misfocusing parameter of step 2.

Ray methods are widely used in depth migration velocity analysis because of its flexibility and low cost compared to wave-equation methods. However, in areas where the highfrequency approximation is violated and the wavefields become too complex, the dependence of velocity-model definition on Kirchhoff migration might result in inaccurate velocity models. This may lead to unreliable results even if wave-equation migration is used to produce the final image. In this case, it is therefore desirable to use wave-equation methods during depth migration velocity analysis.

Importance

Oil exploration is being conducted in geologically complex areas. In the last ten years, huge oil reserves have been continuously discovered below salt bodies in Gulf of Mexico. Recently, it was announced 5 billion barrels of light oil reserve in sub-salt offshore Santos Basin, Brazil. These sub-salt oil discoveries can be partially credited to depth imaging.

As a reliable image in depth is strongly dependent on a accurate velocity model, velocity estimation is one of the most important problems in the exploration community. Velocity model estimation by wave-equation methods is a very active area of research. To be applied in a production environment, however, its cost is still a limiting factor. Therefore, one of the challenges exploration geophysics faces is to decrease the cost of wave-equation methods, specially wave-equation migration, for velocity analysis.

Proposal

Data reduction has a great impact on the final cost of migration. For wave-equation migration velocity analysis, one can take advantage of the linearity of the wave propagation to linearly combine shots using phase encoding functions (Romero et al., 2000; Liu et al., 2006). Decomposition and combination of data in plane-waves is capable of decreasing the amount of data by one order of magnitude in 2D problems (Liu et al., 2006).

Aiming at a significantly higher data reduction, Biondi (2006, 2007) presented an extension of the exploding reflector model to generate prestack data suitable for wave-equation migration velocity analysis. This method uses the exploding reflector concept (Loewenthal et al., 1976) to synthesize areal sources and areal receivers along the entire survey at the surface. The initial conditions for the modeling are prestack images obtained with wave-equation migration.

Here, I propose a hybrid approach to wave-equation migration velocity analysis under the framework of the prestack exploding reflector modeling with one-way extrapolator. By using one-way operators, phase encoding methods can be implemented to adequately reduce input data to migration. I will also investigate the application of modified imaging conditions to attenuate crosstalk due to data reduction. I show that by selecting a few reliable reflectors in conjunction with crosstalk attenuation strategies makes prestack exploding reflector modeling applicable to field data. A question to be addressed is how wave-equation velocity inversion performs under this horizon-based prestack exploding reflector model.

This research will contribute to make wave-equation migration velocity estimation of practical use.

INTRODUCTION

In depth-imaging projects, specially in geologically complex areas, the definition of an accurate velocity model is decisive for achieving a reliable final image. Much effort is expended on building a geologically reasonable velocity model. As velocity model building for depthimaging projects is strongly dependent on Kirchhoff algorithms, however, in areas where ray theory fails, inaccurate updates of the velocity model and even divergence of the solution can occur.

Wave-equation methods have the potential to overcome the limitations ray-based methods face when the velocity model is complex with sharp and irregular boundaries. In waveequation methods, velocity model updates are obtained by minimizing a measure of image misfocussing. The differential semblance optimization (DSO) (Symes and Carazzone, 1991; Shen et al., 2003, 2005) and wave-equation migration velocity analysis (WEMVA) are the most promising wave-equation methods. In particular, WEMVA (Sava, 2004; Sava and Biondi, 2004a,b) measures the misfocusing in a non-automated way by picking a velocity parameter that best flattens angle-Domain common image gathers (ADCIGs) after residual prestack depth migration or, alternatively, residual depth moveout. In the following I give a brief review of WEMVA.

Wave-equation migration velocity analysis

Sava (2004) and Sava and Biondi (2004a,b) introduced WEMVA which iteratively updates the velocity model using an image perturbation, 1I, as a function of the slowness parameter, ρ D \frac{s}{s_{0}}. This parameter is obtained by interactively picking the values that best flatten reflectors on residual migrated ADCIGs. s and s_{0} are the trial slowness and the reference slowness, respectively. The slowness perturbation, 1S D s − s_{0}, is determined by minimizing the objective function:

J (1S) D ||1I − L1S ||^{2}.  (1)

In equation 1, the image perturbation is obtained by applying a linearized version of the residual prestack depth migration onto the reference image, I_{0}, using the interpreted slowness parameter, ρ. Doing so, one ensures that the Born approximation is not violated when the velocity changes are too large. The operator L is the chaining of an imaging operator and the Born scattering operator, S, represented by:

S(U_{0},1S) ≈ −i \frac{dk_{z}}{ds}|_{sDs_{0}}1S1zU_{0}.  (2)

Note that the operator S depends on the reference wavefield ,U_{0}, and on the slowness perturbation, 1S. The term \frac{dk_{z}}{ds}|_{sDs_{0}} is given by:

\frac{dk_{z}}{ds}|_{sDs_{0}} D \frac{ω^{2}s_{0}}{\sqrt{ω^{2}s_{0}^{2}− |k|^{2}}} ,  (3)

where k is the horizontal wavenumber and ω is the angular frequency.

The inversion of the local information given by the image perturbation into global slowness updates is performed by an iterative algorithm, such as conjugate gradients. The image perturbation is back-propagated into slowness perturbation along wave-paths.

In 3D, WEMVA has considerable storage and computational costs (Sava, 2004). Every iteration involves roughly the cost of two wave-equation migrations. The practical implementation of wave-equation migration velocity analysis, therefore, is still a challenge. There are opportunities for improvements like automatic picking of the velocity parameter (Biondi, 2008) and reduction of the storage and computational cost by decreasing the amount of input data to migrate.

Previous data reduction schemes

After the seminal work of Schultz and Claerbout (1978), in which they introduced the concept of synthesis of plane waves, different strategies have been formulated to reduce the computational cost of wave-equation migration by combining shot records based on the linearity of wavefield propagation. The combination of shots drastically reduces the cost of wave-equation migration but generates crosstalk between uncorrelated shots during imaging. The crosstalk can be attenuated by using sufficient linear combinations of shots. In the case of plane-wave migration, if the set of synthesized plane-waves adequately samples the range of propagation angles present in the original data the crosstalk is largely attenuated (Liu et al., 2006).

Rietveld et al. (1992) used the idea of predefined wavefields of arbitrary shape at the subsurface from Schultz and Claerbout (1978) to synthesize areal sources at the surface. These areal sources convolved with shots synthesize areal shots to be migrated using an areal shot migration scheme. Rietveld and Berkhout (1992) further developed the method to prestack migration using a target-oriented controlled illumination strategy. The phase encoding migration of Romero et al. (2000) linearly combine shots using phase functions different from the plane-wave. They use random phase encoding which disperses the crosstalk energy throughout the image and linear phase encoding which shifts the crosstalk outside the migration domain. They show that the migration results deteriorate as the number of encoded shots increases.

The combination of more than one phase encoding technique achieves additional reduction of the computational effort. Sun et al. (2002) and Zhang et al. (2003) show good results on Marmousi data using mixed phase encoding and sign-opposite phase encoding on previously plane-wave encoded data. Liu et al. (2006) provide a general framework to evaluate planewave decomposition in prestack source plane-wave migration and show the equivalence of shot profile and plane-wave migrations. Additionally, the use of plane-wave migration enables the use of different coordinate systems suited to image steep dips (Shan and Biondi, 2004; Shan et al., 2007).

These previous works synthesize a smaller dataset by combining the recorded shot records. Under a different perspective, Biondi (2006, 2007) introduces the prestack exploding reflector modeling. This method synthesizes source and receiver wavefields along the entire survey at the surface, in the form of areal data, starting from a prestack migrated image cube represented by subsurface-offset domain common-image gathers (SODCIGs). For the case of migration velocity analysis, the aim is to generate a small dataset by combining several individual modeling experiments, while maintaining the kinematics necessary to perform migration velocity analysis. As in the other data reduction schemes, crosstalk is also generated during imaging. I will present strategies to attenuate it.

Summary

In the next section, I will first present the prestack exploding reflector model and discuss the crosstalk problem and attenuation methods. Particularly, I show that in the presence of several reflectors the combination of random phase encoding and a depth-window around selected reflectors virtually eliminates the crosstalk. Finally, I will point out the future work as well as the time-line.

WORK COMPLETED

Prestack exploding reflector modeling

The basic idea behind the prestack exploding reflector model (Biondi, 2006, 2007) is that, under the single scattering assumption, the synthesis of 2D source and receiver wavefields at the surface can be described by:

\begin{matrix}S(x,ω) D G(z_{ξ},x_{ξ} − h_{ξ};x,z D 0,ω) ∗ I_{s}(z_{ξ},x_{ξ},h_{ξ}), \\ R(x,ω) D G(z_{ξ},x_{ξ} C h_{ξ};x,z D 0,ω) ∗ I_{r}(z_{ξ},x_{ξ},h_{ξ}), \end{matrix}  (4)

where S(x,ω) is the source wavefield at z D 0, R(x,ω) is the receiver wavefield at z D 0, I_{s}(z_{ξ},x_{ξ},h_{ξ}) is the prestack image used as initial condition to model the source wavefield, I_{r}(z_{ξ},x_{ξ},h_{ξ}) is the prestack image used as initial condition to model the receiver wavefield, and G(z_{ξ},x_{ξ} ± h_{ξ};x,z D 0,ω) represents an operator which extrapolates the wavefields from the subsurface to the surface. Here, I use one-way extrapolators. Additionally, h_{ξ} is the subsurface offset, z_{ξ} is depth, ω is the temporal frequency, x_{ξ} is the spatial position of the selected prestack image, and x is the spatial coordinate in the data space coinciding with x_{ξ}.

The prestack images used as initial condition for the source and receiver wavefield extrapolation are supposed to be dip-independent gathers computed by changing the dip along the offset direction according to the apparent geological dip (Biondi, 2007). In Appendix A, this transformation is described in detail.

To illustrate how the prestack exploding reflector modeling works, I consider a simple 0.75km deep flat reflector in a 2km/s constant velocity model. Figure 1 shows the result of shot-profile migration using a too slow velocity of 1.8km/s. The front panel represents the zero-subsurface offset section and the panel on the right is a SODCIS taken at 2km. Notice, on the SODCIG panel, that the energy does not focus at zero subsurface offset and the reflector is mapped at a shallower depth.

Figure 1: Shot-profile migration using a slower velocity. nodipshtpr [CR]

In Figure 2a is shown the areal source and Figure 2b the areal receiver resulting from the modeling of a single SODCIG as initial condition located at 1.75km. Since the source wavefield is backward propagated, it occurs at negative times. The modeling was performed with the same slower velocity used in migration.

The modeled areal data is to be used in a subsequent areal shot migration using an updated velocity model. However, in this example I show areal shot migration results with the same slower velocity, just for clarification. Figure 3a presents the SODCIG taken at 1.75km after the areal shot migration. For comparison, in Figure 3b is shown the SODCIG resulting from the shot profile migration.

Figures 4a and 4b show the angle-domain common-image gather (ADCIG) computed from data in Figures 3a and 3b, respectively. Notice that the kinematic information relative to migration with too low velocity was maintained.

Figure 2: Areal data modeled from a single SODCIG as initial condition: a) areal source; and b) areal receivers. migar1 [CR]

Figure 3: SODCIG from areal shot migration of the modeled data (a) and from shot profile migration of the original data(b). migscig1 [CR]

Figure 4: ADCIG from areal shot migration of the modeled data (a) and from shot profile migration of the original data(b). Notice that the kinematic information was maintained. migacig1 [CR]

To decrease the amount of input data to migration, areal shots need to be combined into super areal shots after the upward propagation. This can be achieved by regularly selecting individual experiments and adding them up into their respective super areal shots, after being upward propagated. Pairs of super areal source and super areal receiver are to be used as the source and receiver wavefields, respectively, in areal shot migration.

To illustrate the potential data reduction the prestack exploding reflector model can achieve in 3D, I use the geometry of the 3D field data Valenciano (2008) imaged with shot-profile migration:

• number of shots: 10,000 (200 along inline and 50 along crossline)

• migration aperture: 140 ilines and 330 crosslines

• model space: 50 ilines and 150 crosslines

• number of subsurface offsets: 31

Consider generating 3D areal data at every 10th SODCIGs in the inline and crossline directions. So, the number of super areal shots would be 100. As will be shown later, more realizations improve the quality of the areal shot migration when using random phase encoding during the modeling. Consider 4 realizations, so the final number of super areal shots would be 400.

The cost of modeling one areal shot is compared to the cost of the areal migration without the application of imaging condition. Therefore, the total cost to obtain an areal shot migrated image (including the modeling of the areal data) should roughly correspond to the migration of 800 super areal shots. As the aperture of the prestack exploding reflector modeling and areal shot migration spans the whole migrated output space, the relative cost comparison with shot-profile migration is difficult. The cost comparison is easier to be done with plane-wave migration as the migrated domain is the same. For instance, Shan (2008) used 2,700 planewaves in a 3D plane-wave migration.

In a migration velocity estimation workflow, the areal data is likely to be modeled once. Therefore, in this example, the cost of subsequent migrations should correspond to 1/7 of plane-wave migration.

Sava (2004) pointed out that WEMVA has a huge storage cost because the wavefields need to be saved on disk. As in the prestack exploding reflector model the number of areal shots is about to be 1 or 2 orders of magnitude less than that of the original dataset, the storage costs are expected to be dramatically decreased.

As I show next, due to the data reduction, crosstalk arises from correlation of unrelated events during imaging. Crosstalk degrades the image, obscuring the kinematic information we are interested in to perform migration velocity analysis.

Crosstalk attenuation

Wavefield propagation is a linear process. This allows to combine wavefields to decrease the size of input data to migration. However, during imaging, unrelated source and receiver wavefields crosscorrelate, generating crosstalk. Conventionally, to attenuate crosstalk some type of phase encoding technique is applied (Romero et al., 2000; Sun et al., 2002; Zhang et al., 2003; Liu et al., 2006).

In the prestack exploding reflector modeling, let us consider the migration of a super areal data with a single pair of individual experiments. The prestack image is formed by crosscorrelating the source, S̃_{m}, and receiver, R̃_{m}, wavefields according to:

Ĩ_{m}(z_{ξ},x_{ξ},h_{ξ}) D ∑_{ω} S̃_{m}^{⋆}(z_{ξ},x_{ξ} − h_{ξ},ω)R̃_{m}(z_{ξ},x_{ξ} C h_{ξ},ω),  (5)

where ^{⋆} represents complex conjugation. If S̃_{m}(x,ω) and R̃_{m}(x,ω) are comprised by two summed areal shots, the image Ĩ_{m}(z_{ξ},x_{ξ},h_{ξ}) will be given by:

Ĩ_{m}(z_{ξ},x_{ξ},h_{ξ}) D I_{1}(z_{ξ},x_{ξ},h_{ξ}) C I_{2}(z_{ξ},x_{ξ},h_{ξ})C

∑_{ω} S_{1}^{⋆}(z_{ξ},x_{ξ} − h_{ξ},ω)R_{2}(z_{ξ},x_{ξ} C h_{ξ},ω) C

∑_{ω} S_{2}^{⋆}(z_{ξ},x_{ξ} − h_{ξ},ω)R_{1}(z_{ξ},x_{ξ} C h_{ξ},ω).  (6)

In equation 6, the last two summation terms represent the crosstalk.

To illustrate the crosstalk problem in the prestack exploding reflector strategy let us make use of a simple 2km/s constant velocity model with two intersecting reflectors, one horizontal and the other dipping 15^{∘}. The dataset is comprised of 200 split-spread shots spaced every 20 m with a maximum offset of 2000 m. Figure 5 shows the shot-profile migration results using the correct velocity. The trace spacing in the migrated result is 20 m and the number of subsurface offsets in the prestack image is 41, ranging from -400m to 400m. The front panel corresponds to the zero-subsurface offset section, the side panel is a SODCIG selected at x=1.5 km. This prestack image is used to model the areal data.

Figure 5: Prestack cube obtained with shot profile migration with the correct velocity. mig08 [CR]

Figure 6 shows super areal sources (top row) and super areal receivers (bottom row) computed with different distances between SODCIGs. Since the areal source wavefield is backward propagated it exists for negative times. The super areal source labeled a) and super areal receiver labeled d) are modeled from SODCIGs separated every 51 traces (1000 m); the super areal source labeled b) and super areal receiver labeled e) are modeled from SODCIGs separated every 41 traces (800 m); and the super areal source labeled c) and super areal receiver labeled f) are modeled from SODCIGs separated every 11 traces (200 m). These distances in number of traces defines the number of super areal shots to be migrated and, consequently, the migration efficiency. As I show next, this also influences the intensity of the crosstalk.

Figure 6: Super areal data computed with different distances between SODCIGs. Top row: source wavefields. Bottom row: receiver wavefields. Left column: areal data computed with a SODCIG distance of 51 traces (a and d). Central column: areal data computed with a SODCIG distance of 41 traces (b and e). Right column: areal data computed with a SODCIG distance of 11 traces (c and f). ardata [CR]

Figure 7 shows the areal shot migration using cross-correlation imaging condition of the super areal data of Figure 6. The number of subsurface offsets is 41. Figure 7a shows the migration of 51 super areal data from Figure 6a and 6d), Figure 7b shows the migration of 41 super areal shots (Figure 6b and 6e) and Figure 7c shows the migration of 11 super areal shots (Figure 6c and 6f).

Biondi (2006) shows that using a decorrelation distance greater than twice the maximum subsurface offset, the crosstalks can be eliminated. This explains why migration of 51 super areal shots shows no crosstalk in the SODCIG. The migration of the data with 41 super areal shots shows an acceptable crosstalk in the SODCIG. However, the crosstalk is strong in the SODCIG of the data with 11 super areal shots. Notice that for all of the three results a strong crosstalk occurs at zero-subsurface offset (x,z-domain).

Figure 7: Areal shot migration using conventional cross-correlation imaging condition of 51 super areal shots (a), 41 super areal shots (b) and 11 super areal shots (c). armigW [CR]

The crosstalk described in the previous paragraph have two distinct origins. Crosstalk in the zero-subsurface offset section results from the cross-correlation of events in the source wavefield with that in the receiver wavefields modeled by the same SODCIG, which herein I call type-1. Crosstalk in non-zero subsurface offsets is related to cross-correlation of events in the source wavefields with that in the receiver wavefields modeled by different SODCIGs, which herein I call type-2.

Biondi (2007) and Guerra and Biondi (2008) show that type-1 crosstalk is strongly attenuated by windowing the wavefields around zero propagation time. This method yields good results if the migration velocity is accurate. In this case, reflectors ‘explode’ at zero time (Claerbout, 1985). However, the window size must be chosen with care if the migration velocity is inaccurate to include the lags where the wavefields correlate.

Here, I describe the use of random phase encoding during the modeling as a way to attenuate crosstalk. This strategy can be combined with the selection of specific reflectors to be modeled for the case of more complex geology. This procedure will be explained at the end of this section.

Crosstalk attenuation with random phase encoding

Ideally, the areal shots pertaining to the same super areal shot should be uncorrelated. A possible way to decrease the correlation between areal shots is by encoding them with a random phase function in such a way that the crosstalk is dispersed throughout the image. Romero et al. (2000) linearly combine shots, using as encoding function a(x_{ξ},ω) D e^{if (x_{ξ},ω)} . f (x_{ξ},ω) is a random sequence varying between (−π,π).

The modeling with f (x_{ξ},ω)-phase encoding synthesizes data according to

\begin{matrix}S_{a}(x,ω) D a(x_{ξ},ω)G(z_{ξ},x_{ξ} − h_{ξ};x,z D 0,ω) ∗ I_{s}(z_{ξ},x_{ξ},h_{ξ}), \\ R_{a}(x,ω) D a(x_{ξ},ω)G(z_{ξ},x_{ξ} C h_{ξ};x,z D 0,ω) ∗ I_{r}(z_{ξ},x_{ξ},h_{ξ}), \end{matrix}  (7)

and the super areal shots of random phase encoded areal sources and receivers, S̃_{n}^{a}(x,ω) and R̃_{n}^{a}(x,ω), respectively, are given by:

S̃_{n}^{a}(x,ω) D ∑_{nD1}^{k} ∑_{iDn,k,N}S_{i}^{a}(x,ω),

R̃_{n}^{a}(x,ω) D ∑_{nD1}^{k} ∑_{iDn,k,N}R_{i}^{a}(x,ω),  (8)

Figures 8a and 8b show (x_{ξ},ω)-random phase encoded super areal source and receiver wavefields, respectively.

Figure 8: Areal source (a) and receiver (b) wavefields computed with random phase encoding, a(x_{ξ},ω). arrphw [CR]

The image, Ĩ_{m}(z_{ξ},x_{ξ},h_{ξ}), after areal shot migration of super areal data comprised by two randomly phase encoded areal sources and receivers is given by:

Ĩ_{m}(z_{ξ},x_{ξ},h_{ξ}) D I_{1}(z_{ξ},x_{ξ},h_{ξ}) C I_{2}(z_{ξ},x_{ξ},h_{ξ})C

∑_{ω} a_{1}^{⋆}a_{2}S_{m}^{a_{1}^{⋆}}(z_{ξ},x_{ξ} − h_{ξ},ω)R_{m}^{a_{2}}(z_{ξ},x_{ξ} C h_{ξ},ω) C

∑_{ω} a_{2}^{⋆}a_{1}S_{m}^{a_{2}^{⋆}}(z_{ξ},x_{ξ} − h_{ξ},ω)R_{m}^{a_{1}}(z_{ξ},x_{ξ} C h_{ξ},ω).  (9)

The cross-correlation between the phase encoding functions a_{1}(x_{ξ},ω) and a_{2}(x_{ξ},ω) in the summation terms must provide random values that disperse the crosstalk throughout the image.

Figure 9 shows the areal shot migration of 11 super areal shots of (x_{ξ},ω)-randomly encoded areal data. The reduction of the type-2 crosstalk in the SODCIG is significant. However, the type-1 crosstalk in the zero-subsurface offset is much less affected. This crosstalk originates from uncorrelated reflections modeled from the same SODCIG as described before. A possible way to attenuate this crosstalk is to decorrelate the exploding reflectors by using a different random function a(z_{ξ},x_{ξ},ω) D e^{if (z_{ξ},x_{ξ},ω)} .

The modeling with a(z_{ξ},x_{ξ},ω)-random phase encoding synthesizes data according to

\begin{matrix}S_{a}(x,ω) D εa(z_{ξ},x_{ξ},ω)G(z_{ξ},x_{ξ} − h_{ξ};x,z D 0,ω) ∗ I_{s}(z_{ξ},x_{ξ},h_{ξ}), \\ R_{a}(x,ω) D εa(z_{ξ},x_{ξ},ω)G(z_{ξ},x_{ξ} C h_{ξ};x,z D 0,ω) ∗ I_{r}(z_{ξ},x_{ξ},h_{ξ}). \end{matrix}  (10)

In equation 10, the ε term is used to control the randomness of the wavefields. As Romero et al. (2000) pointed out, the encoding results deteriorate as the amount of encoded data increases. By letting ε D 1, the migration results become too noisy. A typical value for ε is 0.1.

Figure 10 show the areal shot migration results of 11 super areal shots of randomly encoded areal data using (z_{ξ},x_{ξ},ω)-randomly encoded areal data. The comparison with Figure 9 and Figure 7c shows that encoding in (z_{ξ},x_{ξ},ω) yields a significant reduction of the type-1 crosstalk in the zero-subsurface offset.

The dispersed noise is not negligible. Using more realizations of (z_{ξ},x_{ξ},ω)-randomly phase encoded data, the noise amplitudes can be largely attenuated. Figure 11 shows the superior result of migrating 4 different realizations of 11 (z_{ξ},x_{ξ},ω)-randomly encoded data. Besides virtually eliminating the crosstalk, the speckled noise has much lower amplitudes than in Figure 10.

In the previous examples, I used simple models consisting of two reflectors to illustrate crosstalk attenuation methods. The crosstalk originated by the correlation of reflections modeled from the same SODCIG are fairly well attenuated by the random phase encoding in (z_{ξ},x_{ξ},ω). However, in the presence of several reflectors this type of crosstalk is severe and using only random phase encoding techniques do not perform successfully. So, I introduce the horizon-based prestack exploding reflector model which is the basis for the hybrid approach of wave-equation velocity inversion. I illustrate the horizon-based prestack exploding reflector model using the sediment part of Sigsbee2a (Paffenholz et al., 2002) dataset.

Figure 9: Areal shot migration of randomly encoded data. 11 super areal shots of (x_{ξ},ω)-randomly encoded areal data were migrated. Compare with Figure 7c. arm11phw [CR]

Figure 10: Areal shot migration of randomly encoded data. 11 super areal shots of (z_{ξ},x_{ξ},ω)–randomly encoded areal data were migrated. Compare with Figure 9. arm11phz [CR]

Figure 11: Areal shot migration of randomly encoded data. 4 groups of 11 super areal shots of (z_{ξ},x_{ξ},ω)– randomly encoded areal data were migrated. Compare with Figures 9 and 10. arm11phz4 [CR] Prestack exploding reflector modeling - Sigsbee2a

I use part of the Sigsbee2a dataset (Paffenholz et al., 2002) to show how the crosstalk attenuation methods described previously perform when more reflectors are present. I selected the sediment part to facilitate the analysis. Again, the correct migration velocity was used in all wavefield extrapolations. The areal shots were generated using every 15^{th} SODCIG as initial condition for the modeling. Therefore, 15 super areal shots were input to migration. Here, I show the results of areal shot migration using the conventional imaging condition in frequency and migration of (z_{ξ},x_{ξ},ω)-randomly encoded areal data.

Figure 12 shows the shot profile migration of the sediment portion of Sigsbee2a data. The main feature in the zero-subsurface offset section is represented by a faulted syncline with reflectors dipping in different directions and a fault on the left side. Figure 13 shows the areal shot migration using the frequency-domain imaging condition of the non-encoded data. This image is dominated by type-1 and type-2 crosstalk represented by the extra ‘reflectors’ when comparing with Figure 12, making the true reflectors unrecognizable.

Figure 12: Shot profile migration of the sediment portion of Sigsbee2a dataset. sigshp [CR]

Figure 13: Areal shot migration of the non-encoded data in the sediment portion of Sigsbee2a using frequency-domain imaging condition. Compare with Figure 12. signicW [CR]

I show the areal shot migration of 4 realizations of both (x_{ξ},ω)-random encoding and (z_{ξ},x_{ξ},ω)-random encoding in Figures 14 and 15, respectively. In general, type-2 crosstalk is reasonably attenuated in the subsurface-offset gathers. Type-1 crosstalk in the zero-subsurface-offset section, however, is still severe even in the migration of (z_{ξ},x_{ξ},ω)-randomly encoded Figure 14: Areal shot migration of 4 realizations of (x_{ξ},ω)-randomly encoded data in the sediment portion of Sigsbee2a. signphw4 [CR]

Figure 15: Areal shot migration of 4 realizations of (z_{ξ},x_{ξ},ω)-randomly encoded data in the sediment portion of Sigsbee2a. signphz4 [CR]

data. This is somewhat related to the degradation of the image quality when the number of encoded shots increases, as observed by Romero et al. (2000).

To circumvent the crosstalk problem in Sigsbee2a example, I use a sparse representation of the initial conditions by selecting only a few reliable, high amplitude reflectors in the original SODCIGs to model. To justify this sparse representation of the migrated image, remind that for velocity model updating we can use some key reflectors, generally high signal-to-noise ratio (S/R), to extract the residual-moveout information (Stork, 1992; Kosloff et al., 1996; Marfurt and Duquet, 1999). Herein I call this procedure horizon-based prestack exploding reflector model.

Figure 16 shows two horizons I choose to constrain the modeling. As the imaged reflectors are characterized by a wavelet stretched to depth, I use a small window around the horizon to select the reflector amplitude to input to modeling.

Figure 17 shows the areal shot migration of 4 realizations of (z_{ξ},x_{ξ},ω)-randomly encoded areal data constrained by the horizons. The crosstalk is almost completely eliminated.

Before proceeding to the next section, where I point out the work to be done, I give a rationale for using the horizon-based prestack exploding reflector model in wave-equation velocity inversion.

The main application of WEMVA is velocity model update under complex overburden, for instance, irregular salt bodies. To obtain accurate updates, it is crucial that the velocity Figure 16: Two picked horizons used to constrain the prestack exploding reflector modeling. sigpick2 [NR]

Figure 17: Areal shot migration of 4 realizations of (z_{ξ},x_{ξ},ω)-randomly encoded areal data using horizons. sigphz4 [CR]

model for the overburden is reasonably defined, such that observed defocusing effects are due to velocity errors in the underlying sediments.

The velocity contrast between salt bodies and sediments is commonly strong. But in general, seismic velocity in sedimentary layers is horizontally smooth. This is because sediment deposition is mainly controlled by the gravitational force (McLane, 1999), so the main velocity variations occur in the vertical direction. Strong lateral velocity variations can occur across big geological faults which put rocks of different ages in lateral contact. Even so, velocity varies smoothly within the faulted blocks.

Apart from the geological reasons, velocity uncertainties are generally big because of complex overburden causing uneven illumination of the subsurface, limited acquisition geometry and presence of noise. These factors make the velocity inversion problem ill-posed and the use of regularization, penalizing roughness of the velocity model, is mandatory (Clapp, 2001). Additionally, velocity resolution is limited by the Fresnel zone size (Thore and Juliard, 1999). For these reasons, velocity interpreters generally prefer smooth velocity models in the sediment portion (Guo and Fagin, 2002). This characteristic is frequently observed in most of the depth migration velocity models from Gulf of Mexico.

Theoretically, by using a few reflectors, the slowness perturbations are supposed to be smoother than that resulting from the inversion of the complete image perturbation. However, given that velocity models are generally smooth below the complex overburden, I expect this possible additional smoothness will not be a problem. In the other hand, these few reliable reflectors can possibly give rise to more reliable velocity updates. These are questions I expect to address in this research.

WORK TO BE DONE

In the section "Work Completed", I explained the prestack exploding reflector model and discussed the crosstalk problem. Finally, I introduced the horizon-based approach to prestack exploding reflector model to make the method applicable to field data. The results on a realistic synthetic seismic data confirms that the method is promising.

To further illustrate the good results of the horizon-based prestack exploding reflector model, Figure 18 shows the data from Figure 17 transformed to the angle domain. For comparison, the ADCIG of the original shot profile migration is shown in Figure 19. Notice that the kinematic information is maintained.

Figure 18: Angle-domain data computed from the areal shot migration of 4 realizations of (z_{ξ},x_{ξ},ω)-randomly encoded areal data constrained by the horizons. sigphz4ang [CR]

Figure 19: Angle-domain data computed from the original shot-profile migration. sigshtang [CR]

As the main goal of this research is to explore WEMVA using the prestack exploding reflector concept, in the following I outline the general structure of WEMVA under the framework of the horizon-based prestack exploding reflector model by exploring the topics below.

Some important questions have to be addressed, like:

• WEMVA in the presence of noise,

• WEMVA in the presence of illumination problems.

Other important issue is how crosstalk attenuation performs when using:

• different phase encoding schemes,

• alternative imaging conditions.

These questions can be properly addressed using 2D algorithms. The 3D extension of the algorithms will enable to tackle more interesting problems from the exploration point of view. So, other developments to be made are to extend to 3D:

• the rotation of SODCIGs,

• the prestack exploding reflector modeling algorithm,

• the areal shot migration.

As I showed, the horizon based approach to the prestack exploding reflector model plays an important role in crosstalk attenuation. So, to incorporate this approach in the 3D workflow it will be necessary to:

• investigate 3D residual prestack migration of the picked horizons,

• investigate 3D and 5D automatic horizon picking strategies.

WEMVA in the presence of noise

Similarly to DSO, one potential problem WEMVA may face is the presence of coherent noise. Mulder and ten Kroode (2002) showed that velocity inversion with DSO is very sensitive to the presence of multiples and converted waves. To obtain reasonable results, a 2D-wavenumber filter was incorporated into the DSO objective function.

Another issue is how WEMVA will perform in the presence of the randomness of the wavefields when using random phase encoding. Or, how the scattered wavefield computed by the interaction of the reference wavefields with velocity perturbations is affected by the randomness of the wavefields? Will the image perturbations, resulting from the migration of the scattered wavefield, be stable enough to provide reliable velocity updates?

To deal with the noise problem, a possibility I want to explore in this research is to apply WEMVA on a modeled dataset with only on a few select reflectors with higher S/R.

WEMVA in the presence of illumination problems

Reflectors do not focus around zero-subsurface offset under uneven illumination. When the illumination problem is too severe, shadow zones occur and migration artifacts dominate the migration result (Valenciano, 2008). In the presence of illumination problems, reflectors in ADCIGs can show departures from flatness not related to velocity errors.

One solution to uneven illumination is to pose imaging as an inverse problem (Clapp, 2005; Valenciano, 2008). Because of ill-posedness, this inverse problem needs regularization. However, regularization in the prestack inversion enforces reflectors to be focused at zerosubsurface offset or, alternatively, flatness in the ADCIG. When the velocity is not accurate this regularization can destroy the kinematic information necessary to velocity updates.

Sava (2004) shows that the angular coverage in the subsurface influences the slowness perturbation. The example of WEMVA applied on Sigsbee clearly shows that the quality of the velocity update deteriorates below the salt body.

I will investigate the potential of WEMVA constrained to some specific horizons. This characterizes a hybrid approach to velocity analysis using wave-equation methods. Ideally, these horizons must limit layers with significant differences in the acoustic properties and correspond to reflectors with good S/R.

The horizon-based technique is used to conveniently constrain the null space of reflectiontomography problems by parameterizing the model space with layers (Stork, 1992; Kosloff et al., 1996; Marfurt and Duquet, 1999). Additionally, this representation agrees with the formulation of the velocity determination problem in a target oriented way, considering that a reasonable velocity model is available for shallower depths. Consequently, for velocity update purposes, just the migrated image along these key horizons is sufficient.

Using different phase encoding schemes

One drawback of applying random phase encoding is that it makes the reflectors randomly “explode” at times different from the zero time of the wavefield propagation. Consequently, the time-windowed imaging condition can not be applied.

The linear phase encoding (Romero et al., 2000) aims to shift the crosstalk out of the imaging domain. The crosstalk starts occurring at propagation times greater than the time-shift applied to the wavefields, enabling us to use the time-windowed imaging condition. Another phase encoding schemes to be tested are the mixed-phase encoding and the sign-opposite phase encoding (Sun et al., 2002; Zhang et al., 2003).

Using alternative imaging conditions

Type-2 crosstalk is similar to the one described by Sava (2007). Sava (2007) uses an imaging condition that cross-correlates decomposed source and receiver wavefields as a function of an a priori known local reflector slope at every position and time. If the reflector slope is not available, it is necessary to loop over a range of possible dip angles in a vicinity of the expected reflector dip. This procedure yields good results, but at the expense of a high computational cost. The horizon-based approach can yield the dip angle information to implement this imaging condition in a more efficient way.

The research on phase encoding schemes and alternative imaging conditions is extremely interconnected and I expect to determine the most appropriate combination between them.

3D extension of rotation, modeling and migration algorithms

The 3D extension of the prestack exploding reflector modeling and areal shot migration algorithms along with code parallelization is, in principle, trivial. However, more understanding is needed to formulate the 3D SODCIG rotation in the case of inaccurate migration velocity. That is because, except for the case of common-azimuth migration (Biondi and Palacharla, 1996), when the migration velocity is inaccurate source and receiver rays do not necessarily cross and the co-planarity condition is not respected.

Picking strategies

As already mentioned, using a horizon-based approach to prestack exploding reflector model virtually eliminates crosstalk. Its efficiency, however, may depend on the S/R of the reflectors and on the vertical distance between them. There seems to be, therefore, a vertical decorrelation distance between reflectors. Its determination is a work still to be done.

For the 2D case of a prestack image obtained with an accurate velocity, just two-dimensional information is needed. If the 2D prestack images are computed with an inaccurate velocity model, however, it is necessary to expand the dimensionality of the picks, using the subsurfaceoffset dimension to correctly capture the kinematic deviations from the zero-subsurface offset. In 3D, the picks are five-dimensional. To implement the horizon-based approach in 3D, automatic picking strategies will be needed, at least to pick the original depth migrated data.

Horizon picking after velocity updates are also needed. I intend to use 3D residual prestack depth migration (Sava, 2003) to transform the original 5D picks into updated picks after velocity updates.

3D residual prestack depth migration accounts for remapping, in the wavenumber domain, the vertical wavenumber axis, k_{z}_{0}, of the migrated image into the vertical wavenumber axis, k_{z}, of the residual migrated image according to:

k_{z} D \sqrt{\frac{1}{ρ^{2}} \frac{(k_{z_{0}}^{2} C |k_{m}|^{2})(k_{z_{0}}^{2} C |k_{h}|^{2})}{4k_{z_{0}}^{2}} − \frac{1}{4} (k_{m} C k_{h})(k_{m} C k_{h})}

\sqrt{\frac{1}{ρ^{2}} \frac{(k_{z_{0}}^{2} C |k_{m}|^{2})(k_{z_{0}}^{2} C |k_{h}|^{2})}{4k_{z_{0}}^{2}} − \frac{1}{4} (k_{m} − k_{h}) (k_{m} − k_{h})}  (11)

where k_{m} and k_{h} are the midpoint and subsurface-offset wavenumber vectors, respectively, and ρ D \frac{v_{0}}{v} is the ratio between the initial velocity, v_{0}, and the updated velocity, v.

TOOLS AND DATA

Software

In the proposed method an initial migration must be available. To generate the initial images I plan to use the one-way wave equation shot-profile migration, coded by Brad Artman. This algorithm has been already used in three PhD thesis to migrate 3D data (Artman, 2007; Alvarez, 2007; Valenciano, 2008). This algorithm is therefore considered quite stable and reliable. To perform WEMVA I will run the codes Paul Sava used in his thesis.

2D prestack exploding reflector modeling and 2D areal shot migration algorithms are already tested. Their 3D extension is still a work to be done along with code parallelization.

SEPlib and SEP3D programs complete the suite of programs needed for my research.

Hardware

Even with the data reduction achievable by using the prestack exploding reflector model, 3D modeling and migration algorithms are computationally intensive. Moreover, 3D WEMVA is a big challenge because of the storage and computational requirements.

The computational resources for this research will be provided by SEP’s computer power. Currently, SEP has 1 shared memory machine with 16 nodes and 32Gb RAM (pompei), 3 shared memory machines with 8 nodes and 8Gb RAM (mad, sad and glad), 1 Linux PC cluster with 32 nodes and 2GB RAM (sep200’s) and 1 Linux PC cluster with 32 nodes and 4/8GB RAM (sep400’s). Two additional Linux PC clusters are expected to be configured soon.

I expect to count on this computational resources to develop my research.

Seismic Data

To prove the concepts and test the programs, besides of modeling specific datasets, SEP data library has interesting data to be used.

I will use Sigsbee2a to test 2D algorithms and 3D SEG-EAGE salt model for the 3D extensions.

Petrobras has elected a challenging offshore 3D data acquired over an area of the Santos Basin, Brazil, to benchmark depth-imaging algorithms and velocity-model estimation strategies. The data was depth migrated with a Kirchhoff algorithm and the velocity was estimated using reflection tomography. In spite of the effort the contractor expended to produce an accurate velocity model, Petrobras evaluates that there are opportunities for improvements. Especially, imaging of steep salt-flanks and characterization of a possible allochtonous salt body in the middle of the section along with sediments below.

Along with the seismic data itself, the depth migration velocity model used to migrate data from Figure 20, observer’s report and 9 time-migrated horizons are available. This data could to be used in this research.

Figure 20: Section of the depth migrated data from Santos Basin, Brazil. fig01 [NR]

If this data proves not be adequate for the purposes of this research, because of lack of sediments below the salt, Elf dataset in the North sea also presents challenges to velocity definition. It is also a well known data set and it has been extensively used by SEP students.

TIME-LINE

1. \underline{Summer 2007-2008}

Further investigation on crosstalk attenuation and imaging conditions 2D tests of WEMVA using horizon-based prestack exploding reflector model Investigation on picking strategies 3D migration of SEG-EAGE salt model Course requirements

2. \underline{Autumn 2008-2009}

Further investigation on crosstalk attenuation and imaging conditions Investigation on picking strategies 3D extension of the areal shot migration 3D extension of the prestack exploding reflector model Test 3D algorithms SEG-EAGE 3D salt model Course requirements

3. \underline{Winter 2008-2009}

3D extension of the areal shot migration 3D extension of the prestack exploding reflector model Investigation on picking strategies Test 3D algorithms SEG-EAGE 3D salt model Course requirements

4. \underline{Spring 2008-2009}

3D migration of field data Picking of the field data migration Course requirements

5. \underline{Summer 2008-2009}

3D prestack exploding reflector modeling of field data 3D velocity model update of field data 3D areal shot migration of field data Start writing thesis

6. \underline{Autumn 2009-2010}

3D prestack exploding reflector modeling of field data 3D velocity model update of field data 3D areal shot migration of field data Writing thesis

7. \underline{Winter 2009-2010}

Finalize deliverables for graduation

APPENDIX A

Rotation of Common Image Gathers

In shot-profile migration, SODCIGs are computed by applying the multi-offset imaging condition (Rickett and Sava, 2002), in which the wavefields are shifted horizontally prior to correlation. If the migration velocity is inaccurate, the image point does not coincide with the point where source and receiver rays intersect (Biondi and Symes, 2004). Consequently, in the presence of dip, source- and receiver-ray end points do not lie on the imaged reflector. To correct for the effects these imperfections have on the prestack exploding reflector modeling results Biondi (2007) proposes to apply a preprocessing step before modeling.

In the preprocessing step, a rotation is applied to the CIGs, changing the dips in the z − h plane (depth-offset) based on the apparent dips in the z − x plane (depth−x coordinate). By rotating the CIGs, the upward propagation of source and receiver wavefields starts along the apparent source and receiver ray local directions. The rotation can be implemented in the wavenumber domain (k_{z}− k_{x} − k_{h}).

As Figure A-1 shows, local directions of the source and receiver rays are given by θ_{s} D α −γ and θ_{g} D α Cγ, respectively, where α is the apparent dip and γ is the apparent reflection angle. The aperture angle is related to the dips along the offset direction by (Sava and Fomel, 2000):

tanγ D −\frac{k_{h}_{ξ}}{k_{z}_{ξ}},  (A-1)

where k_{z}_{ξ} is the vertical wavenumber and k_{h}_{ξ} is the offset wavenumber. The apparent geological dip is computed as

tanα D −\frac{k_{x}_{ξ}}{k_{z}_{ξ}},  (A-2)

where k_{x}_{ξ} is the midpoint wavenumber. Using equations A-1 and A-2 and trigonometric expressions, tanθ_{s} and tanθ_{g} are computed as

tanθ_{s}(k_{z}_{ξ},k_{x}_{ξ},k_{h}_{ξ}) D k_{z}_{ξ} \frac{k_{h}_{ξ} − k_{x}_{ξ}}{k_{z}_{ξ}^{2}C k_{x}_{ξ}k_{h}_{ξ}},  (A-3)

tanθ_{g}(k_{z}_{ξ},k_{x}_{ξ},k_{h}_{ξ}) D −k_{z}_{ξ} \frac{k_{h}_{ξ} C k_{x}_{ξ}}{k_{z}_{ξ}^{2}− k_{x}_{ξ}k_{h}_{ξ}} ,  (A-4)

The rotation of the reflectors in the CIGs is realized via a simple remapping of subsurfaceoffset wavenumber according to

k̂_{h}_{ξ}^{s} D −k_{z}_{ξ} tanθ_{s}  (A-5)
k̂_{h}_{ξ}^{g} D k_{z}_{ξ} tanθ_{g}  (A-6)

where k̂_{h}_{ξ}^{s} and k̂_{h}_{ξ}^{g} are the wavenumbers corresponding to the new subsurface-offset axes ĥ_{ξ}^{s} and ĥ_{ξ}^{g}, respectively.

Figure A-1: Geometric representation of the rotation of an image event from horizontal subsurface offset to geologic-dip subsurface offset (from Biondi (2007)). cig-rot [NR]

To illustrate the importance of the SODCIG rotation, I use a simple model of a 25^{∘}-dipping reflector. Data were migrated with a 10 percent slower velocity. The SODCIG at 2km was used as non-rotated (Figure A-2) initial condition and after rotation (Figure A-3) to generate areal data. For the non-rotated case, the initial condition is the same for the source and receiver wavefield modeling.

Figure A-2: Initial condition to model the source wavefield and the receiver wavefield for the non-rotated case. dipshtpr1 [CR]

The modeled areal data were migrated using the correct velocity (2km/s). After migrating with the correct velocity, the energy is expected to focus in the up-dip direction, since it was initially under-migrated (Thore et al., 2002). Figures A-4 and A-5 show the areal shot migration using areal data modeled from non-rotated SODCIG and areal data modeled from rotated SODCIG, respectively. Areal migration using data modeled from rotated SODCIG correctly positions the event and focus the reflection around the zero-subsurface offset as can be seen in the right panel. Areal migration using data modeled from non-rotated SODCIG does not laterally shift the reflector and the reflection energy shows some curvature in the subsurface-offset gather.

Figure A-3: SODCIG after rotation according to the geological dip used as initial condition to model the source wavefield (a) and the receiver wavefield (b). datarot [CR]

Figure A-4: Areal shot migration using areal data modeled from nonrotated SODCIG of Figure A-2. dipnorot [CR]

Figure A-5: Areal shot migration using areal data modeled from rotated SODCIG of Figure A-3. diprot [CR]

REFERENCES

Alvarez, G., 2007, Attenuation of multiples in image space: Ph.D. thesis, Stanford University. Artman, B., 2007, Passive seismic imaging: Ph.D. thesis, Stanford University. Biondi, B. and G. Palacharla, 1996, 3-d prestack migration of common-azimuth data: Geo-

physics, 61, no. 6, 1822–1832. Biondi, B. and W. W. Symes, 2004, Angle-domain common-image gathers for migration ve-

locity analysis by wavefield-continuation imaging: Geophysics, 69, no. 5, 1283–1298. Biondi, B., 2006, Prestack exploding-reflectors modeling for migration velocity analyis: SEP–

124, 45–60. Biondi, B., 2007, Prestack modeling of image events for migration velocity analysis: SEP-131,

pages 101–117. Biondi, B., 2008, Automatic wave-equation migration velocity analysis: SEP-134, pages 65–

78. Claerbout, J., 1985, Imaging the earth’s interior: Blackwell Scientific Publications,Inc. Clapp, R., 2001, Geologically constrained migration velocity analysis: Ph.D. thesis, Stanford

University. Clapp, M. L., 2005, Imaging under salt: illumination compensation by regularized inversion:

Ph.D. thesis, Stanford University. Guerra, C. and B. Biondi, 2008, Prestack exploding reflector modeling: The crosstalk prob-

lem: SEP-134, pages 79–91. Guo, N. and S. Fagin, 2002, Becoming effective velocity-model builders and depth imagers,

part 2—the basics of velocity-model building, examples and discussions: The Leading Edge, 21, no. 12, 1210–1216.

Kosloff, D., J. Sherwood, et al., 1996, Velocity and interface depth determination by tomogra-

phy of depth migrated gathers: Geophysics, 61, no. 5, 1511–1523. Liu, F., D. W. Hanson, et al., 2006, Toward a unified analysis for source plane-wave migration:

Geophysics, 71, 129–139. Loewenthal, D., L. Lu, R. Roberson, and J. Sherwood, 1976, The wave equation applied to

migration: Geophysical Prospecting, 24, 380–399. Marfurt, K. and B. Duquet, 1999, Mapping prestack depth-migrated coherent signal and noise

events back to the original time gathers using fermat’s principle: Geophysics, 64, 934–941. McLane, M., 1999, Sedimentology: Oxford University Press.

Mulder, W. A. and A. P. E. ten Kroode, 2002, Automatic velocity analysis by differential semblance optimization: Geophysics, 67, no. 4, 1184–1191.

Paffenholz, J., J. Stefani, B. McLain, and K. Bishop, 2002, SIGSBEE 2A Synthetic subsalt dataset – Image quality as function of migration algorithm and velocity model error: 64th Annual International Meeting, Eur. Assoc. Geosc. Eng., Extended Abstracts, page B019.

Rickett, J. E. and P. C. Sava, 2002, Offset and angle-domain common image-point gathers for shot-profile migration: Geophysics, 67, 883–889.

Rietveld, W. E. A. and A. J. Berkhout, 1992, Prestack depth migration by means of controlled illumination: Geophysics, 59, 801–809.

Rietveld, W. E. A., A. J. Berkhout, and C. P. A. Wapenaar, 1992, Optimum seismic illumination of hydrocarbon reservoirs: Geophysics, 57, 1334–1345.

Romero, L., D. Ghiglia, C. Ober, and S. Morton, 2000, Phase encoding of shot records in prestack migration: Geophysics, 65, 426–436.

Sava, P. and B. Biondi, 2004a, Wave-equation migration velocity analysis. I. Theory: Geophysical Prospecting, 52, 593–606.

Sava, P. and B. Biondi, 2004b, Wave-equation migration velocity analysis. II. Subsalt imaging examples: Geophysical Prospecting, 52, 607–623.

Sava, P. and S. Fomel, 2000, Angle-gathers by Fourier Transform: SEP–103, 119–130.

Sava, P. C., 2003, Prestack residual migration in the frequency domain: Geophysics, 68, no. 2, 634–640.

Sava, P., 2004, Migration and velocity analysis by wavefield extrapolation: Ph.D. thesis, Stanford University.

Sava, P., 2007, Stereographic imaging condition for wave-equation migration: Geophysics, 72, no. 6, A87–A91.

Schultz, P. S. and J. Claerbout, 1978, Velocity estimation and downward continuation by wavefront synthesis: Geophysics, 43, 691 714.

Shan, G. and B. Biondi, 2004, Imaging overturned waves by plane-wave migration in tilted coordinates: SEG Technical Program Expanded Abstracts, 23, no. 1, 969–972.

Shan, G., R. Clappand, and B. Biondi, 2007, 3d plane-wave migration in tilted coordinates: SEG Technical Program Expanded Abstracts, 26, no. 1, 2190–2194.

Shan, G., 2008, Imaging of steep reflectors in anisotropic media by wavefield extrapolation: Ph.D. thesis, Stanford University.

Shen, P., W. W. Symes, and C. C. Stolk, 2003, Differential semblance velocity analysis by wave-equation migration: SEG Technical Program Expanded Abstracts, 22, no. 1, 2132– Shen, P., W. W. Symes, et al., 2005, Differential semblance velocity analysis via shot profile

migration: SEG Technical Program Expanded Abstracts, 24, no. 1, 2249–2252.

Stork, C., 1992, Reflection tomography in the postmigrated domain: Geophysics, 57, 680– 692.

Sun, P., S. Zhang, and F. Liu, 2002, Prestack migration of areal shot records with phase encoding: 72nd Ann. Internat. Mtg, Soc. Expl. Geophys., Expanded Abstracts, pages 1172–1175.

Symes, W. W. and J. J. Carazzone, 1991, Velocity inversion by differential semblance optimization: Geophysics, 56, no. 5, 654–663.

Thore, P. D. and C. Juliard, 1999, Fresnel zone effect on seismic velocity resolution: Geophysics, 64, no. 2, 593–603.

Thore, P., A. Shtuka, et al., 2002, Structural uncertainties: Determination, management, and applications: Geophysics, 67, no. 3, 840–852.

Valenciano, A., 2008, Imaging by wave-equation inversion: Ph.D. thesis, Stanford University.

Zhang, S., P. Sun, N. Yi, and S. Li, 2003, Prestack migration of areal shot records with mix phase encoding: SEG Technical Program Expanded Abstracts, 22, no. 1, 1031–1034.