Pre-Stack Exploding-Reflector Model

This chapter introduces the pre-stack exploding-reflector model (PERM). PERM uses exploding reflectors as the initial conditions to synthesize data, in a manner similar to the well-known exploding-reflector model (ERM). However, PERM also considers reflectivity as a function of subsurface-offsets, as opposed to the zero-subsurface offset initial condition used by ERM. In this sense, PERM is a generalization of ERM since wave-equation migration of PERM data can generate a pre-stack image, which is not achievable with ERM. The modeling of PERM data can be performed using any wave propagation scheme; here I use the one-way wave-equation. As in the explodingreflector model, the modeling experiments can be combined, potentially decreasing data size in orders of magnitude when compared to the original shot records. Data size reduction is very appealing for migration velocity analysis, especially using wavefield extrapolation methods in 3D projects. After describing the exploding-reflector concept, I generalize it by describing the theory of PERM. The usefulness of PERM to migration velocity analysis, is demonstrated through migration examples.

INTRODUCTION

Migration is applied to seismic data to generate an image of the subsurface. For many years, migration was applied only in the post-stack domain, using the idea of exploding reflectors (??). The simple but powerful concept of the exploding-reflector model (ERM) states that a zero-offset time section can be considered as the recording at the surface of wavefields generated by simultaneous explosions of all points in the subsurface. The strength of the explosions is proportional to the reflection coefficient and, because the wavefield propagates from the subsurface to the surface, to obtain correct kinematics the medium velocity must be halved.

Migration in the post-stack domain assumes that data is transformed to some approximation of the zero-offset acquisition geometry. Because of the required transformation to zero-offset, many algorithms have been developed like dip moveout (DMO) (???), migration to zero-offset (MZO) (??), and common-reflection surface (??). DMO and MZO can be considered very particular cases of the more general azimuth moveout (AMO) (?), which must be applied to 3D data prior to commonazimuth migration (?).

In fairly simple geology, ERM perfectly matches the kinematics of the zero-offset data that would have been recorded with coincident source and receivers at the surface. This means that for ERM to hold, it is necessary that the downgoing path from the source to a point in the subsurface be the same as the upgoing path from the point in the subsurface back to the receiver. However, this assumption is often invalid in areas of geological complexity. Post-stack imaging below salt bodies, for instance, can give rise to numerous migration artifacts, in the presence of prismatic and strongly refracted waves (?). Therefore, post-stack migration does not produce reliable images in the presence of strong lateral velocity variation.

In areas of complex geology pre-stack depth migration becomes mandatory not only for imaging purposes but also for velocity estimation. In such areas, migration by wavefield extrapolation has been widely used to produce the final image because it properly handles complex distortions of the wavefields. However, due to the high computational cost, wavefield extrapolation methods are rarely used to estimate the migration velocity model in 3D projects (?), where ray-based methods are the industry standards. In addition to the lower cost, ray-based methods are very flexible with respect to strategies for defining the velocity model (????). But despite their advantages, ray-methods do not satisfactorily describe complex wave propagation in the presence of large lateral velocity contrasts. In this case, a more complete description of the wavefield complexities is needed, and therefore we face the challenge of decreasing the cost of migration velocity analysis by wavefield extrapolation while maintaining its robustness.

A typical way of decreasing the cost of wavefield extrapolation is to reduce the amount of input data. Data size can be reduced by selecting a smaller number of shots. Figure ?? shows pre-stack images computed with the multi-offset imaging condition (?) using 380 shots 24 m apart (Figure ??a) and 23 shots 384 m apart (Figure ??b). The shot records are modeled with the one-way wave equation, using a smoothed version of the Marmousi velocity model. The same velocity model is used to migrate the data. The right panel is the zero-subsurface-offset section and the left panel is the subsurface-offset-domain common-image gather (SODCIG) selected at x = 2760m on the right panel. Notice the slanted lines in the SODCIG of Figure ??b representing the angles according to which reflectors are illuminated. When using all the shots, the slanted lines constructively interfere and the energy is focused at zero-subsurface-offset (Figure ??a). perm01width=a) Shot-profile migration of 380 off-end shots 24 m apart. b) Shot-profile migration of 23 off-end shots 384 m apart. Both images were computed with the correct velocity model. Notice slanted lines present in Figure ??b.

By selecting only a few shots, the migration output shows poor angular coverage. An alternate way of reducing data size is by combining shot-profiles into areal shots based on the linearity of wavefield propagation. Data is multiplied by a comb function and stacked to originate one areal shot. The comb function is shifted until all shots are selected. The original angular coverage is maintained if the period of the sampling function is sufficiently big (Figure ??a). However, if many shots are combined, crosstalk is generated since unrelated shots and receiver wavefields are cross-correlated during imaging (Figure ??b). Compare Figures ??a and ??b with Figure ??a. Crosstalk can completely overwhelm the reflectors, and the kinematic information for migration velocity analysis can be lost. perm02width=a) Shot-profile migration of 95 areal shots resulting from the combination of 4 shot-profiles 2256 m apart. b) Shot-profile migration of 23 areal shots resulting from the combination of 16 shot-profiles 564 m apart. Both images were computed with the correct velocity model. Notice crosstalk occurring periodically in the SODCIG of Figure ??b.

The combination of wavefields is implicit in ERM. When reflectors are allowed to simultaneously explode, linearity of wavefield propagation is evoked to combine wavefields initiated at every point in the subsurface. However, crosstalk is not generated in post-stack migration based on ERM since no cross-correlation of wavefields is performed because the imaging condition is a simple summation over frequency, which extracts the image at zero time of wavefield propagation (?). Therefore, it is natural to consider the combination of wavefields in any generalization of ERM.

Combination of wavefields is exploited by the prestack-exploding reflector model (PERM) (?) to significantly reduce the data size used in migration velocity updates. Like ERM, PERM uses the concept of exploding-reflectors to propagate wavefields. However, instead of considering reflectivity only as a function of the spatial coordinates, reflectivity is also parameterized as a function of subsurface-offsets. Its elementary idea is to synthesize data necessary to correctly image a single reflector of an isolated SODCIG, so that migration of PERM data shows correct kinematics needed for performing velocity updates.

An interesting feature of PERM is that it can be used in a target-oriented way, since the modeled data contains all the information necessary to image a predefined region of the subsurface. This concept is present in different methods, such as controlled illumination (??) and common-focus point (CFP) technology (???).

In this chapter, I introduce and further develop PERM. It will be shown that data synthesized with PERM has the kinematic information necessary to perform migration velocity analysis. Although, in this thesis, migration velocity analysis is performed using wave-equation tomography, ray-based tomography can also be used in conjunction with PERM data. Before introducing the theory of PERM, I will first describe the exploding-reflector concept of which PERM is a generalization. The usefulness of PERM for migration velocity analysis will be illustrated by comparing areal-shot migration of PERM wavefields with that of shot-profile migration for simple bidimensional models. The extension to a tridimensional (3D) medium is straightforward. In 3D, it will be shown that under the common-azimuth approximation (?), PERM drastically decreases data size.

EXPLODING-REFLECTOR MODEL

To describe the exploding-reflector model (ERM) let us analyze the modeling of seismic data under the acoustic approximation. Here, we consider the Born or singlescattering approximation. This consideration leads us to a linear operator whose adjoint is the migration operator.

Let us start with the constant-density acoustic wave equation for a single temporal frequency ω

(∇^{2} + ω^{2}s^{2}(x)) P (x, ω) = 0,  (1)

where ∇^{2} is the Laplacian operator, s(x) is the slowness field, and P (x, ω) is the wavefield. Note that x = (x, y, z) is the vector of spatial coordinates. By introducing reflectivity as r(x) = 1− \frac{s^{2}(x)}{s_{0}^{2}(x)}, where s_{0}(x) is a smooth background slowness, equation 1 can be written as

(∇^{2} + ω^{2}s_{0}^{2}(x, ω)) P (x, ω) ≈ ω^{2}s_{0}^{2}(x)r(x)P (x, ω).  (2)

The total wavefield P (x, ω) can be considered as the sum of an incident background wavefield P_{0}(x, ω) with a perturbed or scattered wavefield ΔP (x, ω). Therefore, we can write

(∇^{2} + ω^{2}s_{0}^{2}(x)) ΔP (x, ω) ≈ ω^{2}s_{0}^{2}(x)r(x)P (x, ω),  (3)

given that the background wavefield P_{0}(x, ω) is the solution of equation 2 using the background slowness s_{0}(x). Notice that equation 3 is a non-linear relationship between r(x) and ΔP (x, ω). A linear operator can be derived by using Green’s function and the Born (weak scattering) approximation. The background Green’s function G_{0}(x′, x, ω) is the solution of equation 1 using the background slowness in the presence of a point source δ(x − x′) at x′ = (x′, y′, z′):

\{ \begin{matrix}(∇^{2} + ω^{2}s_{0}^{2}(x)) G_{0}(x′, x, ω) = 0 \\ G_{0}(x′, x = x′, ω) = δ(x − x′) \end{matrix} .  (4)

Multiplying equation 4 by [ω^{2}s_{0}^{2}(x)r(x)P (x, ω)], integrating with respect to x′ over the volume V in the subsurface, and comparing the result with equation 3, we see that

ΔP (x, ω) ≈ _{∫V} ω^{2}s_{0}^{2}(x′)r(x′)G_{0}(x′, x, ω)P (x′, ω)dx′.  (5)

Under the weak scattering assumption, the total wavefield P (x′, ω) can be approximated by the background wavefield P_{0}(x′, ω), and the linear relationship between the scattered wavefield ΔP (x, ω) and reflectivity r(x) reads

ΔP (x, ω) ≈ _{∫V} ω^{2}s_{0}^{2}(x′)r(x′)G_{0}(x′, x, ω)P_{0}(x′, ω)dx′.  (6)

Using equation 6 and assuming that the background slowness field and the reflectivity distribution are known, the scattered wavefields ΔP (x_{s}, x_{r}, ω) measured by receivers at x_{r} = (x_{r}, y_{r}, z_{r}) due to sources located at x_{s} = (x_{s}, y_{s}, z_{s}) are

ΔP (x_{s}, x_{r}, ω) ≈ _{∫V} ω^{2}s_{0}^{2}(x)G_{0}(x_{s}, x, ω)r(x)G_{0}(x, x_{r}, ω)dx,  (7)

which amounts to convolving the source Green’s function with the reflectivity, then with receiver Green’s function, and summing the contributions of all points within V. By using one-way propagators, we downward propagate the source wavefield, convolve it with the reflectivity, and upward propagate this result up to the surface.

Zero-offset data is obtained by selecting from ΔP (x_{s}, x_{r}, ω) traces with x_{s} = x_{r},

ΔP (x_{s}, ω)|_{x_{s}=x_{r}} ≈ ∫ ω^{2}s_{0}^{2}(x)r(x) G_{0}(x_{s}, x, ω)|_{x_{s}=x_{r}} G_{0}(x, x_{s}, ω)dx.  (8)

ERM synthesizes zero-offset data by initiating virtual sources located on the reflectors. The fundamental consideration is that the downgoing and upgoing rays of a zero-offset source and receiver pair follow the same path. The exploding-reflector wavefield is propagated with half of the actual background velocity and the virtual sources explode with strength proportional to the reflection coefficient. Assuming that the background slowness and the reflectivity distribution are known, the wavefield generated by the exploding-reflector model P_{ERM}(x, ω) propagates according to the following one-way wave-equation:

\{ \begin{matrix}(\frac{∂}{∂z} + i\sqrt{ω^{2}s^{2}(x) − |k|^{2}}) P_{ERM}(x, ω) = r(x) \\ P_{ERM}(x, y, z = z_{max}, ω) = 0 \end{matrix} ,  (9)

ERM is kinematically correct only if source and receiver Green’s functions are equal for x_{s} = x_{r} at a point x in the subsurface. Notice that this condition is unlikely to be fulfilled in the presence of strong lateral velocity variations (?). Moreover, ERM does not model multiples and prismatic reflections (?).

In equation 9, the reflectivity r(x) acts as the initial condition for the wavefield propagation. Notice that a migrated image can replace the reflectivity in the initial condition. Taking into consideration the imaging principle (?) and the computation of pre-stack images by wave-extrapolation methods (??), ERM, as initially formulated, implicitly assumes that all the seismic energy is perfectly focused at the zerosubsurface offset and that it is sufficient to parameterize the migrated image only as a function of the position vector x. This implies that the slowness field is accurate and illumination of the subsurface is sufficiently good. When this is not the case, the migrated image must be described not only as a function of x but also as a function of the subsurface offset or aperture angle.

The generalization of ERM proposed in this thesis does not aim at computing prestack data. Rather, it synthesizes data whose migration with wavefield-extrapolation methods enables the computation of a pre-stack image. PERM uses the ERM concept with a poorly focused pre-stack image as initial condition to propagate wavefields. The modeled data is, potentially, orders of magnitude smaller than the original shot records and contains all necessary kinematic information to update the velocity model using ray-based tomography or, as used in this thesis, wave-equation-based tomography. This characterizes the pre-stack-exploding-reflector model, which will be described in the next section.

PRE-STACK-EXPLODING-REFLECTOR MODEL

The fundamental idea of PERM is to model data which describes the correct kinematics of an isolated SODCIG. Many shot records contribute to form the image at a point in the subsurface. Therefore, to model data using the conventional one-way modeling, we would perform several modeling experiments consisting of downward continuing the source wavefield initiated by point sources at the surface, convolving the propagated source wavefield with the SODCIG, and upward continuing the convolution result up to the surface. As we do not know beforehand which shots contribute to forming the image at a point in the subsurface, we would have to model every shot originally present in the original dataset.

Ideally, instead of performing many modeling experiments, we would like to synthesize a small amount of data with the condition that migration has the same kinematics as the initial SODCIG. This can be achieved by using a strategy similar to ERM. To compute SODCIGs using the new modeled data, we need to synthesize source and receiver wavefields.

The modeling of source and receiver wavefields can be carried out by any wavefieldcontinuation scheme. Here, we use the following one-way wave equations:

\{ \begin{matrix}(\frac{∂}{∂z} − i\sqrt{ω^{2}s_{0}^{2}(x) − |k|^{2}}) D_{P}(x, ω; x_{m}) = I_{D}(x_{m}, h) \\ D_{P}(x, y, z = z_{max}, ω; x_{m}) = 0 \end{matrix} ,  (10)

and

\{ \begin{matrix}(\frac{∂}{∂z} + i\sqrt{ω^{2}s_{0}^{2}(x) − |k|^{2}}) U_{P}(x, ω; x_{m}) = I_{U}(x_{m}, h) \\ U_{P}(x, y, z = z_{max}, ω; x_{m}) = 0 \end{matrix} ,  (11)

where I_{D}(x_{m}, h) and I_{U}(x_{m}, h) is the isolated SODCIG at the horizontal location x_{m} for a single reflector, suitable for the initial conditions for the source and receiver wavefields, respectively. The subsurface-offset h can be parameterized as h = (h_{x}, h_{y}, h_{z}). In this thesis, when describing 2D problems h = (h_{x}), and for the 3D case, h = (h_{x}, h_{y}) without considering the computation of the vertical subsurface offset h_{z} as introduced by ?. The initial conditions are obtained by rotating the original unfocused SODCIGs according to the apparent geological dip of the reflector. For dipping reflectors, this rotation maintains the velocity information needed for migration velocity analysis. The rotation of SODCIGs will be described later in this chapter. If the initial condition has energy focused at the zero-subsurface offset or if no pre-stack information is available, the pre-stack image can be parameterized only by its spatial coordinates, and PERM is equivalent to ERM.

Let us now illustrate the generation of PERM data synthesized from a single SODCIG. We start with a pre-stack image computed with shot-profile migration of 401 split-spread shots at every 10 m and maximum offset of 2250 m, using a 10% lower velocity (Figure ??). The model consists of one reflector at 750 m depth embedded in a medium with constant velocity of 1000 m/s. The pre-stack image has 81 subsurfaceoffsets ranging from -400 m to 400 m. Notice the poor focusing of energy around the zero-subsurface offset due to inaccurate velocity.

The SODCIG at x_{m} = 0 m was used as the initial condition for modeling the corresponding pair of PERM source and receiver wavefields using the same inaccurate velocity. The wavefields are upward propagated according to equations 10 and 11. The PERM data is shown in Figure ??. Notice that the receiver wavefield (Figure ??a) occurs at positive times while the source wavefield (Figure ??b) occurs at negative times. According to the imaging principle (?), reflectors explode at time zero. This is the time at which the source wavefield impinges on the reflector. Because the receiver wavefield exists after the source wavefield has reached the reflector, the areal receiver data U_{P}(x, y, z = 0, ω; x_{m}) is upward propagated forward in time. For the same reason, the areal source data D_{P}(x, y, z = 0, ω; x_{m}) is upward propagated backward in time.

refpl01width=Shot-profile migration of 401 split-spread shots 10 m apart with a 10% slower velocity. The model consists of a horizontal reflector embedded in constant velocity of 1000 m/s.

refpl02width=Data synthesized by PERM having as the initial condition the SOD-CIG at x_{m} = 0 m. a) The receiver wavefield. b) The source wavefield.

Areal-shot migration of data from Figure ?? with the same inaccurate velocity shows kinematics at near subsurface-offsets similar to that of the original shot-profile migration (compare Figures ?? and ??). For farther subsurface-offsets, however, energy is not adequately imaged. This can be easily explained by analyzing how the image I(x, h) is formed when applying the multi-offset imaging condition

I(x, h) = ∑_{ω} ∑_{x_{s}} D^{∗}(x − h, x_{s}, ω)U(x + h, x_{s}, ω).  (12)

In equation 12, D(x − h, x_{s}, ω) and U(x + h, x_{s}, ω) are the source and receiver wavefields, respectively, shifted laterally by h. The image is formed by summing the contributions of all shot positions x_{s} and all frequencies. Notice that ‘^{∗}’ stands for the complex conjugate. From equation 12 we can see that farther subsurface-offsets get contributions from the cross-correlation of wavefields within a neighborhood around x. Specifically for the PERM case, farther subsurface-offsets of the areal-shot migrated image at x_{m} (the location of the isolated SODCIG) also get contributions from the cross-correlation of wavefields generated at SODCIGs within a neighborhood around x_{m}. Therefore, to ensure that the areal-shot migrated image has kinematics at all subsurface-offsets similar to that in the original isolated SODCIG, we need to model PERM data from a set of SODCIGs around x_{m}. From equation 12, it is clear that the maximum distance at which the wavefields still contribute to form the image at the maximum subsurface-offset is given by 2(max |h_{x}| , max |h_{y}|).

refpl03width=Areal-shot migration of PERM data shown in Figure ?? with a 10% slower velocity. By comparing with Figure ?? we see that far subsurface-offsets are not properly imaged.

The areal-shot migration of PERM data synthesized by isolated SODCIGs within the interval (−2max |h_{x}| , 2max |h_{x}|) is shown in Figure ??. By using more data, energy is adequately imaged at far subsurface-offsets (compare with Figures ?? and ??). To further understand the behavior of the areal-shot migrated image, let us examine the reflection angle-domain common-image gathers (ADCIGs) (?). ADCIGs computed from the SODCIGs of Figures ??, ?? and ?? attest to the more accurate imaging when migrating data modeled from the set of SODCIGs around x_{m} (Figure ??). Although the ADCIG from the the image computed with a single pair of PERM data (Figure ??b) shows reasonable kinematics, the amplitude of wide-aperture angles is weaker than that of the original ADCIG (Figure ??a). Notice that the amplitude behavior of the ADCIG computed with several areal shots from SODCIGs within the neighborhood of x_{m} (Figure ??c) better matches that of the original isolated SODCIG.

refpl04width=Areal-shot migration of PERM data having a set of isolated SOD-CIGs around x_{m} = 0 m as the initial condition with a 10% slower velocity. By comparing with Figure ?? we see that the kinematics of far subsurface-offsets is properly recovered. refpl05width=ADCIGs selected at x_{m} = 0 m. a) Computed from the shot-profile migration; b) computed from the areal-shot migration of one pair of PERM data modeled from the SODCIG at x_{m} = 0 m; and c) computed from the areal-shot migration of pairs of PERM data modeled from a set of SODCIGs around x_{m} = 0 m. Notice that although the kinematics is similar, the amplitudes in c) better match those of a).

If, instead of using the incorrect migration velocity, we input the correct migration velocity to the areal-shot migration, energy nicely focuses at zero-subsurface offset (Figure ??). This property will be used to perform migration velocity updates in Chapter 3.

refpl06width=Areal-shot migration of PERM data having a set of isolated SOD-CIGs around x_{m} = 0 m as the initial condition with the correct velocity. Energy nicely focuses at zero-subsurface offset.

In the previous examples we saw that PERM data contains all the kinematic information needed to perform migration velocity analysis. However, by carefully examining the migration results of the original data and PERM data we can see that the latter has stronger side lobes.

The stronger side lobes are due to the squaring of the wavelet. It occurs because PERM source and receiver wavefields contain the waveform of the original image. Let us consider, for simplicity, a plane reflector so that the initial conditions are the same for modeling source D_{P} and receiver U_{P} wavefields. For only one frequency, PERM can be described by

D_{P}(ξ) = L_{0}(ξ, x − h)I(x, h),  (13)

and

U_{P}(ξ) = L_{0}(ξ, x + h)I(x, h).  (14)

The operator L is the recursive upward propagation of the wavefields from the subsurface (represented by the the coordinates x) up to the surface (represented by the coordinates ξ). This operator is also a function of the subsurface-offset because of the way the pre-stack image is injected into the modeling. The subscript of L denotes that the wavefield propagation is performed with the background slowness s_{0}(x) used to migrate the original shot-profiles. The source wavefield is upward propagated backwards in time and the receiver wavefield is upward propagated forward in time.

The pre-stack image I_{P} is computed with PERM data, at every depth level, using the multi-offset imaging condition

I_{P}(x, h) = D∗_{P}(x − h)U_{P}(x + h).  (15)

The wavefields in depth are recursively downward propagated according to

D∗_{P}(x) = L′_{1}(ξ, x)D_{P}(ξ),  (16)

and

U_{P}(x) = L′_{1}(ξ, x)U_{P}(ξ).  (17)

where ′ stands for the complex conjugate. Note that in equations 16 and 17 the subscript of the operator L indicates the use of a different migration velocity. The source wavefield is downward propagated forward in time, and the receiver wavefield is downward propagated backwards in time. Inserting equation 13 into equation 16 and equation 14 into equation 17, and substituting the expressions for the wavefields in depth into equation 15 we get

I_{P}(x, h) = {L′_{0}(ξ, x − h)L_{1}(ξ, x − h)I(x, h)}
× {L′_{1}(ξ, x + h)L_{0}(ξ, x + h)I(x, h)} . (18)

From equation 18, we can see that the pre-stack image I_{P} is approximately a squared version of the original image. This means that in addition to the stronger side lobes, reflectors in the PERM image will always have positive polarity.

For the simple case of a horizontal reflector in a constant velocity medium, we have shown that migration of PERM data produces images with the same kinematics as the shot profile migration. Now, we introduce a dipping reflector in the same constant background velocity medium. In the presence of a non-zero geological dip, a pre-processing of the initial conditions is necessary to obtain correct kinematics. This pre-processing step is represented by a rotation of the pre-stack image according to the apparent geological dip.

Dip-independent initial conditions

Shot-profile and areal-shot migrations by wavefield extrapolation computes pre-stack images by means of the multi-offset imaging condition (?), in which source and receiver wavefields are laterally shifted prior to time correlation. However, the shift between wavefields might not be restricted to the horizontal direction. For instance, vertical shifts of the wavefields produce the vertical-subsurface offset gathers which provide reliable velocity information in the presence of steep dips (?).

Ideally, wavefields should be shifted along the geological dip direction. According to ?, SODCIGs computed this way do not suffer from image-point dispersal in the presence of dip and migration velocity errors. The image-point dispersal makes events with different reflection angle from the same reflection point in the subsurface be imaged at different locations.

The image-point dispersal in 2D is illustrated in Figure ?? for the case of migrating with a slower velocity (Figure ??a) and higher velocity (Figure ??b) than the true velocity. Let us consider constant velocity in the vicinity of the image point, so source and receiver rays are straight for simplicity.

When migration velocity is too low, the reflector is imaged at a smaller depth. The image point computed with horizontal shifts of the wavefields I_{hx} is shifted down-dip with respect to the image point computed with shifts along the apparent geological dip I_{hg}. The geological dip is called apparent because of the migration velocity error. The point I is where source and receiver rays cross at an angle which is twice the apparent reflection angle γ. Source and receiver rays cross deeper than the image points, causing events curve down in the SODCIG.

When migration velocity is too high, the reflector is imaged at a greater depth. The image point I_{hx} is shifted up-dip with respect to I_{hg}. Source and receiver rays cross shallower than the image points, causing events curve up in the SODCIG. pdispwidth=Geometry for the computation of SODCIGs. Source, receiver and image points are labeled with S, R and I, respectively. The subscript hx corresponds to subsurface offsets computed with horizontal shift. The subscript hg corresponds to subsurface offsets computed by shifting along the apparent geological dip α. a) too low velocity, and b) too high velocity. Modified from ?.

Generating SODCIGs along the geological-dip direction overcomes the problem of the image-point dispersal. However, it is computationally demanding since wavefields must be stored at various depths and, in addition, accurate dip information is difficult to obtain, especially when events cross because of velocity inaccuracy.

? point out that, at least in the first order, the reflection-angle domain is immune to image-point dispersal. This is because the SODCIG to ADCIG transformation shifts events to the line connecting I and I_{hg} in Figure ?? at the same image point shared by all the reflection angles.

In the presence of dip, to accurately model PERM data it is crucial that the initial conditions are free of image-point dispersal to guarantee that all the energy of a point in the subsurface is contained by the correspondent SODCIG injected to modeling. Since SODCIGs along the geological dip are not easily computed, can we pre-process the SODCIGs computed with horizontal shifts of the wavefields such that they are transformed into a good approximation of the SODCIGs along the geological dip?

To answer this question, let us first examine the angle relationships in Figure ??. The angles γ + α and γ − α are the source and receiver ray angles, respectively. They are the propagation direction of the wavefields locally at the image point. Notice that α and γ are related to slopes in the pre-stack image according to

tan α = −\frac{dz_{m}}{dx_{m}}  (19)

and

tan γ = −\frac{dz}{dh_{x}} ,  (20)

where the subscript m refers to the local nature of the relationship. The solutions of the differential equations 19 and 20 define slant-stack paths, which allow us to transform the 2D pre-stack image I(x, z, h_{x}) into I(x, z, α, γ) by angle decomposition according to the following integrals:

I(x, z, α, γ) = ∫_{x_{m_{i}}}^{x_{m_{f}}} ∫_{−h_{x}}^{h_{x}} W (x_{m} − x)\frac{dI(x, z, h_{x})}{dz}dx_{m}dh_{x}|_{z=z_{h}+h_{x} tan γ}  (21)
z = z_{m} + x_{m} tan α

(22)

where the derivative with respect to z is to recover the correct phase. The local window W (x_{m} − x) is used in the local slant-stack integral on x_{m}, being defined as

\{ \begin{matrix}1, x_{m_{i}} ≤ x ≤ x_{m_{f}} , \\ 0, elsewhere \end{matrix}

where x_{m_{i}} = x_{m} − \frac{x_{w}}{2} and x_{m_{f}} = x_{m} + \frac{x_{w}}{2} , x_{w} is the width of the local window.

Again, using simple trigonometry, we have

tan(γ + α) = \frac{tan γ + tan α}{1 − tan γ tan α},  (23)

tan(γ − α) = \frac{tan γ − tan α}{1 + tan γ tan α}.  (24)

To align the initial conditions with the geological dip, we need to change the dip along the subsurface-offset axis according to the apparent geological dip, yielding the new subsurface offset h̃_{x_{s}} and h̃_{x_{r}} for the initial conditions of the modeling of source and receiver wavefield, respectively. This is accomplished by solving the following differential equations:

tan(γ + α) = −\frac{dz}{dh̃_{x_{s}}} ,  (25)

tan(γ − α) = \frac{dz}{dh̃_{x_{r}}} .  (26)

The solutions of equations 25 and 26 define new slant-stack operations which, in combination with equations 23 and 24, reduces the dimensionality of the decomposed pre-stack image (equation 27) by transforming I(x, z, α, γ) into I_{D}(x, z, h̃_{x_{s}}) and I_{U}(x, z, h̃_{x_{r}}).

In 3D, we must also consider rotation of the cross-line offsets according to the apparent geological dip in the cross-line direction in addition to the in-line rotation. To apply the rotation in 3D, we need to assume that source and receiver rays are coplanar, such that they cross. In 3D, this assumption might not be valid when migrating with an inaccurate velocity in complex overburden. The co-planarity assumption is also present in the formulation of the common-azimuth migration (?), which turned out to provide sufficiently accurate images in areas of complex geology.

Instead of performing a true 3D slant-stack, which can be computer memory intensive, we can take advantage that a 3D slant-stack can be split into two 2D slantstacks on the in-line and cross-line directions. This can be easily seen by using Fourier transforms to conveniently express a generic 3D slant-stack integral

\begin{aligned}U(ϕ_{1}, ϕ_{2}, ϕ_{3}) & = & ∫ ∫ u(z, x, y)dx dy |_{t=ϕ_{1}+ϕ_{2}x+ϕ_{3}y} \\ & = & ∫ ∫ u(ϕ_{1} + ϕ_{2}x + ϕ_{3}y, x, y)dx dy \end{aligned}  (27)

as

U(ϕ_{1}, ϕ_{2}, ϕ_{2}) = ∫ dϕ_{1}e^{−ikϕ_{1}} ∫ dxe^{−iϕ_{2}x} ∫ dye^{−iϕ_{3}y} ∫ dze^{ikz}u(z, x, y).  (28)

where ϕ_{1}, ϕ_{2} and ϕ_{3} are, respectively, the output z, the slanting parameter corresponding to x, and the slanting parameter corresponding to y.

To illustrate the generation of dip-independent initial conditions 801 split-spread shots 10 m apart with maximum offset of 3250 m were modeled with velocity of 1000 m/s and migrated with a 10% slower velocity (Figure ??). The model has a 20^{∘} dipping reflector and a horizontal reflector at a depth of 2500 m. The SODCIG located at 0 m was used as the initial condition for the modeling of PERM data without applying the pre-processing described above. Since PERM models one event of one isolated SODCIG, the dipping reflector and the horizontal reflector originate two different pairs of PERM wavefields. This means that reflectors used in the modeling need to be interpreted in the pre-stack volume. The pair of source and receiver wavefields for the dipping reflector are shown in Figure ??. dip01width=Shot-profile migration of 801 split-spread shots 10 m apart with velocity 10% slower than the true velocity. The model is represented by a 20^{∘} dipping reflector and a horizontal reflector at depth of 2500 m embedded in constant velocity of 1000 m/s. dip02width=Data synthesized by PERM having as the initial condition the dipping reflector in the SODCIG at x_{m} = 0 m. a) The receiver wavefield. b) The source wavefield.

Since the modeling of a single non-rotated SODCIG carries no dip information, migration of the corresponding PERM data using the correct velocity does not shift events laterally as can be seen in Figure ??. As expected, the horizontal reflector focuses at zero-subsurface offset. However, notice how the dipping reflector still presents a residual curvature. Migration of PERM data from SODCIGs within a neighborhood around x_{m} = 0 m is shown in Figure ??. Again, the residual curvature is present and any migration velocity analysis using this result will lead to wrong velocity updates. This residual curvature is a result of not having corrected the image-point dispersal. Unless stated, SODCIGs in the figures are selected at x = 0 m.

dip03width=Areal-shot migration of PERM data shown in Figure ?? using the correct velocity. The horizontal reflector is focused at zero-subsurface offset, but the dipping reflector shows residual curvature. dip04width=Areal-shot migration with correct velocity of PERM data having a set of isolated SODCIGs around x_{m} = 0 m as the initial conditions. As in Figure ??, the horizontal reflector is focused at zero-subsurface offset, but the dipping reflector shows residual curvature.

The rotation was applied to the image in Figure ??, and the new PERM data was modeled using the initial conditions shown in Figure ??. Notice that the initial condition for modeling the source wavefield (Figure ??a) and the initial condition for modeling the receiver wavefield (Figure ??b) have the dipping event oriented in opposite directions in the SODCIG. The rotation changes neither the horizontal reflector nor the-zero subsurface offset, as can be seen on the right panels. dip05width=Initial conditions for modeling a) source and b) receiver wavefields. The dipping reflector is oriented in opposite directions in the SODCIG. Rotation does no affect the horizontal reflector nor the-zero subsurface offset, as can be seen on the right panels.

The source and receiver wavefields for the dipping reflector after rotation are shown in Figure ??. The events in Figures ??a and ??b are shown in the same areal shot for illustration only. Actually, they pertain to different areal shots. dip06width=Dipindependent PERM data for the dipping reflector from the rotated SODCIG at x_{m} = 0 m. a) The receiver wavefield. b) The source wavefield.

Areal-shot migration of dip-independent PERM data is shown in Figure ??. Notice that the segment of the dipping reflector is shifted laterally with respect to that of the horizontal reflector. Since the dip-independent wavefields carry information about the dip of the reflector, the observed reflector movement is now consistent with migration with a higher velocity. dip07width=Areal-shot migration with the correct velocity of dip-independent PERM data having the rotated the SODCIGs at x_{m} = 0 m as the initial condition. Compare with Figure ??. The dipping reflector is now focused, in contrast to the image in Figure ??, where it shows residual curvature.

Migration with the correct velocity of dip-independent PERM data modeled from a set of SODCIGs in a neighborhood around x_{m} = 0 m confirms the correctness of the rotation (Figure ??). The focusing of the dipping reflector around zero-subsurface offset is greatly improved when compared with Figure ??. The corresponding AD-CIGs confirm the more consistent move-out after rotation (Figure ??). Note the residual move-out in the angle gather corresponding to the image computed with wavefields with non-rotated initial conditions (Figure ??a) and how much flatter is the image from wavefields computed with the proposed rotation (Figure ??b). dip08width=Areal-shot migration with correct velocity of dip-independent PERM data having a set of rotated SODCIGs around x_{m} = 0 m as the initial conditions. Compare with Figure ??. The focusing of the dipping reflector is greatly improved when using the rotated initial conditions. dip09width=0.66ADCIGs of images computed with the correct migration velocity using PERM data having: a) non-rotated initial conditions, and b) rotated initial conditions. Note the residual move-out in a) and the flatter response in b).

In the example with two reflectors, we initiated each modeling experiment from one isolated SODCIG of one single reflector to compute a pre-stack image restricted to a certain region in the output space. However, depending on the number of reflectors and the size of the prestack image, this procedure can generate a dataset even larger than the original shot-profiles, defeating the original purpose of PERM, which is to synthesize a smaller dataset to be used in migration velocity analysis. We see next, that combination of modeling experiments can be used to decrease the size of PERM data.

Combination of modeling experiments

In the previous examples, if we were to fully image the reflectors, the number of areal shots might be twice as many as the original shot-profiles. To decrease the number of modeling experiments, linearity of wave propagation can be used to combine isolated SODCIGs and inject them simultaneously into one single model experiment as

\{ \begin{matrix}(\frac{∂}{∂z} − i\sqrt{ω^{2}s_{0}^{2}(x) − |k|^{2}}) D̂_{P}(x, ω; x̂) = Î_{D}(x̂, h) \\ D̂_{P}(x, y, z = z_{max}, ω; x̂) = 0 \end{matrix} ,  (29)

and

\{ \begin{matrix}(\frac{∂}{∂z} + i\sqrt{ω^{2}s_{0}^{2}(x) − |k|^{2}}) Û_{P}(x, ω; x̂) = Î_{U}(x̂, h) \\ Û_{P}(x, y, z = z_{max}, ω; x̂) = 0 \end{matrix} ,  (30)

where Î_{D}(x̂, h) and Î_{U}(x̂, h) are the combination of SODCIGs for a single reflector to be used as the initial conditions for the modeling of combined wavefields, D̂_{P}(x, ω; x̂) and Û_{P}(x, ω; x̂), the source and receiver wavefields, respectively. Notice that x̂ is a series of evenly spaced points on the reflector, where the SODCIGs are selected. In 3D, the selection of SODCIGs can be thought of as the multiplication of the pre-stack image by a 2D comb function. This sampling function is shifted laterally to select new set of SODCIGs to initiate the modeling of another pair of combined wavefields. After shifting along one period of the sampling function in the x and y directions, all the points on the reflector are used in the modeling. Consequently, the number of modeling experiments equals the number of lateral shifts of the sampling function.

The choice of the sampling period determines how much crosstalk will be present in the migrated image, in a manner similar to the shot-profile migration illustrated in the introduction of this chapter. To obtain a crosstalk-free image, the sampling period needs to be set to a distance such that wavefields initiated at different SOD-CIGs do not correlate. As previously shown when discussing how to compute PERM images with kinematics similar to that of the original shot-profile migration, PERM wavefields generated from SODCIGs within the interval 2(max |h_{x}| , max |h_{y}|) still contribute to the image at the central SODCIG. For the same reason, no crosstalk is generated during migration if the period of the sampling function is larger than that interval. Since the focusing of energy in the SODCIG is velocity error dependent, so is the the period of the sampling function. Therefore, for small velocity errors a small sampling period can be used and, consequently, a smaller number of combined modeling experiments is needed.

The combination of modeling experiments is illustrated in Figure ??. We model PERM data starting with the rotated images from the previous section and combine SODCIGs into sets using the sampling period of 41, 81 and 163 SODCIGs. As can be seen in Figures ??a and ??b, crosstalk occurs in the SODCIG according to a period which is half of the sampling period. Recalling that the number of subsurfaceoffsets is 81, it is expected that no crosstalk will occur when migrating the set with sampling period of 163 SODCIGs (Figure ??c). In this case, since each reflector is used separately in the modeling, the total number of areal shots is 326, which is less than half of the original shot profiles. comb01width=0.85Areal-shot migration of PERM data synthesized from sets of SODCIGs selected with sampling period of:

a) 41, b) 81, and c) 163 SODCIGs. Notice that no crosstalk is generated when the sampling period is larger than twice the maximum absolute subsurface-offset.

To further decrease the number of modeling experiments, we can combine the reflectors in the initial condition. By doing so, two events are present both in the source wavefield and in the receiver wavefield, as in Figure ??. During migration, in addition to the expected cross-correlations, cross-correlation of wavefields from different reflectors will occur, generating crosstalk between reflectors. Areal-shot migration of PERM data synthesized from a set of SODCIGs with sampling period of 163 is shown in Figure ??. Both reflectors are simultaneously injected in the initial condition. Notice the presence of reflector crosstalk labeled with ‘C’. comb02width=0.85Areal-shot migration of PERM data synthesized from a set of SODCIGs selected with sampling period of 163. The two reflectors are simultaneously injected to the model. Notice the reflector crosstalk resulting from the cross-correlation of the wavefields from the horizontal reflector with that from the dipping reflector.

To illustrate the reflector crosstalk generation, snapshots of the wavefield propagation are shown in (Figure ??). Both wavefields are represented in the same panel and signaled with ‘D’ (from downgoing) for the source wavefield and ‘U’ (from upgoing) for the receiver wavefield. The panels on the left show the wavefields in the time-depth domain selected at horizontal positions where reflector crosstalk occurs in Figure ??. The panels on the right are taken at the propagation time when the wavefields cross on the left panels. The crossing times vary from 0.14 s to 1.24 s, for the SODCIGs from right to left, respectively. The crossing times are a function of the distance between reflectors and the propagation velocity. In Figure ??a notice that for propagation times less than 0.14 s minus the period of the wavelet in time no crosstalk will occur. comb03width=0.80Snapshots of the wavefield propagation of the Figure ??. The panels on the left are selected at a horizontal position where the crosstalk occurs in Figure ??. The panels on the right are taken at the propagation time when the wavefields cross on the left panel.

Crosstalk can be mitigated if we take advantage of the imaging principle, which states that the reflector image is formed at zero time of wavefield propagation. This is especially true when using migration velocity is accurate. When the migration velocity is inaccurate, the focusing of the image departs from the zero time, and this can potentially be used to update the velocity model (??). As Figure ?? shows, crosstalk is formed at times different from zero. Therefore, if the wavefields are crosscorrelated within a time window centered at zero, reflector crosstalk might not occur depending on the distance in depth between the reflectors (?). The time-windowed imaging condition for a single areal shot reads

I Pw(x, h) = ∑_{\frac{−tw}{2} ≤t≤ \frac{tw}{2}} F^{−1} [D^{∗}(x − h)] F^{−1} [U(x + h)] ,  (31)

where t_{w} is the length of the time window. When using one-way propagators, the wavefields are inverse Fourier transformed to time by F^{−1}. Migration using the timewindowed imaging condition of equation 31 is shown in Figure ??. The length of the time window is 0.06 s which avoids the cross-correlation of events from different reflectors. The image is completely free of reflector crosstalk. comb04width=0.9Areal-shot migration using the time-windowed imaging condition (equation 31). The reflector crosstalk is completely eliminated.

We saw that by carefully combining the modeling experiments the data size is decreased while the correct kinematics is maintained, which is important for migration velocity analysis. However, combining the modeling experiments using a decorrelation distance between events does not achieve a significant data reduction. In this case, data reduction depends on the number of subsurface-offsets which are necessary to capture all the relevant velocity information. In the example using the time-windowed imaging condition, the number of independent experiments was reduced five times if compared with the original dataset. Data reduction techniques like plane-wave decomposition, for instance, could lower that number ten times. As we will see in Chapter ??, further data reduction can be achieved by using the phase-encoding technique (?) to linearly combine the modeling experiments. This will enable us to use a shorter sampling period of SODCIGs, and also injecting more than one reflector in the modeling.

Although PERM theory was developed in 3D, all the examples so far were 2D. Next, it is discussed an 3D example under the common-azimuth approximation (?). It will be shown that in this case the SODCIGs in the y direction can be continuously sampled and the number of modeling experiments will be determined only by sampling period in the x direction, drastically decreasing data size.

3D-PERM from common-azimuth migrated images

In the way PERM is formulated there is no restriction on the dimensionality of the pre-stack image used as the initial condition for the modeling, which means that if the original data have sufficient cross-line offsets as in the acquisition geometries with wide-range of azimuths (???), the initial conditions are a five-dimensional hypercube on x, h_{x} and h_{y}.

To synthesize PERM data starting with the five-dimensional initial conditions such that no crosstalk is generated during migration, the minimum number of modeling experiments is 4n_{h_{x}}n_{h_{y}}, where n_{h_{x}} and n_{h_{y}} are the number of subsurface offsets in x and y directions. Considering usual parameters, the number of modeling experiments would be as low as some hundreds or few thousands. This data reduction is very substantial if we compare, for instance, with data reduction achieved by 3D-plane-wave migration. Using plane waves, to obtain artifact-free SODCIGs due to the lack of illumination from some propagation directions we need to migrate roughly 2000 plane waves. This means that 3D-PERM data size can be one order of magnitude smaller than 3D-plane wave data.

Despite the recent good migration results obtained in geologically complex areas using wide-azimuth data, narrow-azimuth acquisition is still the industry standard. Narrow-azimuth data can be efficiently imaged by common-azimuth wave-equation migration (CAM) (?). That is because the dimensionality of the pre-stack wavefields is decreased by assuming zero-cross-line offset in contrast with the full-azimuth migration. That does not mean that the cross-line offset wavenumber is zero. Rather, its asymptotic approximation is function of the the in-line midpoint and in-line offset wavenumbers. Therefore, instead of a five-dimensional hypercube CAM images are four-dimensional hypercubes on x and h_{x}.

Because of the lower dimensionality of CAM images, when using them as the initial conditions to synthesize PERM data, the SODCIGs in the cross-line direction can be sampled continuously as depicted in Figure ??b in contrast with the fivedimensional initial conditions for the full azimuth case of Figure ??a. Recall that PERM is equivalent to ERM if energy is focused at the zero-subsurface offset as well as if this is the only available subsurface offset.

The continuous sampling of SODCIGs in the cross-line direction yields one more order of magnitude of data reduction. So, under the common-azimuth approximation, 3D-PERM data size can be two orders of magnitude smaller than 3D-plane wave data. cam01width=1.1The initial conditions for synthesizing PERM data from CAM images can be specified as in a) because no pre-stack information exists in the crossline direction, in contrast with the full azimuth situation in b).

To illustrate the validity of the above assumptions, a split-spread data with maximum offset of 1587.5 m was computed using 3D-Born modeling (?) on a 30^{∘} dipping reflector with 45^{∘} azimuth with respect to the acquisition direction, which is aligned with the in-line direction. The number of in-lines and cross-lines is 96 each spaced at every 25 m. The offset interval is 25m. The velocity used in the modeling is the 1D function v(z) = (1500 + 0.5z) m/s.

The Born data was input to CAM with a 5% slower velocity. Migration results can be seen in Figures ??a and ??b for SODCIGs positioned at (x = 750 m,y = 600 m) and (x = 750 m,y = 1000 m), respectively. The panel on the left is the SODCIG which contains 21 subsurface-offsets ranging from −250 m to 250 m. The panel in the middle is the in-line at the zero-subsurface offset, and y = 600 m (Figure ??a) and y = 1000 m (Figure ??b). The panel on the right is the cross-line at the zerosubsurface offset, and x = 750 m. cam02width=0.9Common-azimuth migration of 3D-Born data modeled from a 30^{∘} dipping reflector with 45^{∘} azimuth with respect to the acquisition direction. The panel in the middle is the in-line at the zero-subsurface offset, and y = 600 m (Figure ??a) and y = 1000 m (Figure ??b). The panel on the right is the cross-line at the zero-subsurface offset, and x = 750 m.

In the common-azimuth regime, the computation of the dip-independent initial conditions is performed by simply rotating the SODCIGs in the in-line direction, since no cross-line offset is computed in migration. The underlying assumption that source and receiver rays cross, which makes the rotation valid in 3D, implies that they are coplanar. This assumption is in agreement with the co-planarity condition present in the formulation of CAM.

Source and receiver wavefields were modeled using continuous sampling along the cross-line direction and sampling period of 48 in the in-line direction. This period is sufficient to not generate crosstalk during the areal-shot migration, given that the number of subsurface-offsets of the pre-stack image is 21. The synthesized 3D receiver wavefield is shown in Figure ??. The left panel is the in-line at y = 1200 m, the right panel is the cross-line at x = 1400 m, and the top panel is the time-slice at t = 0.5 s. cam03width=0.93D-PERM receiver wavefield. The left panel is the in-line at y = 1200 m, the right panel is the cross-line at x = 1400 m, and the top panel is the time-slice at t = 0.5 s.

The 3D migration of the 48 areal shots with the 5% slower velocity is shown in Figures ??a and and ??b for SODCIGs positioned at (x = 750 m,y = 600 m) and (x = 750 m,y = 1000 m), respectively. For facilitating the comparison with the CAM images of Figures ??a and ??b, the polarity of the areal-shot migrated image is inverted due to the squaring of the wavelet. The kinematics of the SODCIGs computed with PERM wavefields matches that of the SODCIGs computed with CAM, especially between −200 m and 200 m. This enables the use of 3D PERM wavefields computed from CAM images in optimization of migration velocity, as will be shown in Chapter ?? for a 3D survey from the North-Sea. cam04width=0.93D-areal-shot migration of PERM data. The panel in the middle is the in-line at the zero-subsurface offset, and y = 600 m (Figure ??a) and y = 1000 m (Figure ??b). The panel on the right is the cross-line at the zero-subsurface offset, and x = 750 m. Compare with Figure ??.

CONCLUSIONS

In this chapter we saw that wavefields synthesized by PERM are able to provide migrated images with correct kinematics while decreasing data size. Data reduction is achieved by combining the modeling experiments and is controlled by the number of subsurface offsets which will be computed during areal-shot migration of PERM data. Recall that SODCIGs in the initial conditions must be separated by a distance at least twice the maximum absolute subsurface-offset value to prevent crosstalk. Implicit to PERM is the fact that reflectors must be identified such that reflector crosstalk is avoided during migration. 3D Pre-stack interpretation can be cumbersome, but it allows, for instance, avoiding the use of reflectors with low signal-to-noise ratio in migration velocity estimation. Moreover, in commercial softwares for migration velocity estimation reflector picking is a usual and almost entirely automated procedure.

While in 2D data reduction with PERM is comparable to that of the plane-wave decomposition, in 3D it can be of one order of magnitude less if cross-line subsurfaceoffsets are to be computed. If the initial conditions are computed with commonazimuth migration, PERM data size can potentially be two orders of magnitude less than 3D plane-wave data.