# Focusing analysis

Series Investigations in Geophysics Öz Yilmaz http://dx.doi.org/10.1190/1.9781560801580 ISBN 978-1-56080-094-1 SEG Online Store

The idea of a velocity analysis that is based on differential solutions of the scalar wave equation first was introduced by Doherty and Claerbout. They used the 15-degree finite-difference migration algorithm and worked with single CMP gathers. Gonzalez-Serrano and Claerbout later extended the wave equation velocity analysis to slant-midpoint coordinates and worked with linearly moveout-corrected CMP gathers.

## Method

The method discussed here operates in the Fourier transform domain using the exact form of the double-square-root (DSR) operator. Mathematical details of this method are found in Section E.7. Figure 5.4-23 summarizes the main computational steps involved in this migration velocity analysis based on wavefield extrapolation.

1. Starting with the prestack data in midpoint y, off-set h, and two-way event time t in the unmigrated position, represented by the wavefield P(y, h, τ = 0, t) at the surface τ = 0, perform 3-D Fourier transform. The variable τ is associated with the direction of wave extrapolation, and is related to depth z by τ = 2z/υ, where υ is the medium velocity.
2. Specify an extrapolation velocity function that only varies vertically, υ(τ) and apply the extrapolation operator exp(–iωDSRτ/2) to compute the extrapolated wavefield in the transform domain P(ky,kh,τ,ω) from the surface wavefield in the transform domain P(ky,kh,τ = 0,ω).
3. To obtain the zero-offset image, sum over the offset wavenumber, and thus obtain P(ky,h = 0, τ, ω).
4. Apply 2-D inverse Fourier transform to obtain the zero-offset image P(y, h = 0,τ,t). The image below a midpoint y is contained in the t – τ plane.
5. Perform mapping of the variables as described in Section E.6 from τ to υ. The velocity information is given by the envelope of the velocity volume of data P(y, h =0, τ = t, υ).

We now demonstrate the procedure outlined in Figure 5.4-23 using a synthetic data set. Figure 5.4-24 shows two common-offset sections over a number of point scatterers buried in a constant-velocity earth, where υ = 3000 m/s. Using a constant velocity for extrapolation, υe = 3000 m/s, the t – τ image plane was produced for each midpoint. Two such planes corresponding to midpoints 1 and 5 denoted in Figure 5.4-24 are shown in Figure 5.4-25. The υ–τ planes (Figure 5.4-26) then were generated from the t – τ image planes by the mapping procedure described in Section E.7. Peak amplitudes for all events occur at the correct medium velocity (3000 m/s). We expect the diffraction events in Figure 5.4-23 to migrate to the apexes beneath midpoint 1, where the point scatterers are located. Note that in Figure 5.4-25, almost all the energy is in the image plane corresponding to midpoint 1; just five midpoints away, at midpoint 5, the migrated energy is very low.

How do we interpret the t – τ image planes? If we used the true medium velocity in downward extrapolation, then, according to the imaging principle, we would see all the events along the diagonal τ = t, the image line, on the image plane. This happens in Figure 5.4-25, because a 3000-m/s extrapolation velocity was used, which is just the velocity used in generating the model in Figure 5.4-23. Any displacement of peak energy from the t = τ image line means that the velocity value used for downward extrapolation differs from that of the event. This displacement is also the basis for mapping from the tτ image plane to the υ – τ plane by equation (E-77).

This mapping is investigated further with the modeled data set shown in Figure 5.4-27, in which velocity increases with depth. In Figure 5.4-28b, note that the top and middle events fall to the left of the image line suggesting that the velocity used in extrapolation (υe = 3000 m/s) is greater than the velocities associated with these events. The bottom event falls on the image line, implying that its velocity is nearly the same as the extrapolation velocity. These observations are confirmed in the corresponding υ – τ planes in Figure 5.4-29. While true stacking velocity values for the three events are 2700, 2850, and 3000 m/s, the velocities interpreted from Figure 5.4-29b are about 2500, 2800, and 3000 m/s. Thus, the migration-based velocity estimate for the shallow event is in error by approximately 8 percent.

To determine the reason for the velocity error, we will consider a migration-based velocity analysis of our synthetic data example that does not involve the approximate mapping step. Figure 5.3-30a shows a CMP gather from midpoint 1 in the zero-dip region of the depth-variable velocity model associated with the constant-offset sections in Figure 5.4-27. The migration velocity analysis on this gather (Figure 5.4-30b) was done by extrapolating the surface wavefield P(kh, ω, τ = 0) repeatedly with different constant velocities in steps of Δτ = Δt (the sampling rate). The zero-offset trace from each attempt was collected after this effort, abandoning the rest of the migrated CMP gather.

Interpretation of the velocity analysis in Figure 5.4-30b reveals correct stacking velocities for the three events, including the shallowest. Clearly the error observed in Figure 5.4-29 is attributable to the mapping (equation (E-100). Note that the error does not occur because of depth variability of velocity, but instead, because the single extrapolation velocity used differed from the medium velocity. The conventional velocity analysis for midpoint 1 of this model data set is shown in Figure 5.4-30c for comparison. Note the familiar NMO stretching that is apparent in the shallow event. In other respects, both the results (Figures 5.4-30b and 5.4-30c) are comparable.

The departure of an event on the t – τ image plane from the t = τ image line is measured by the quantity Δτ as depicted in Figure 5.4-31a. In some practical implementations, the t – τ image plane is mapped onto the plane of Δτ versus τ as depicted in Figure 5.4-31c to determine the rms velocity υ(τ) for time migration from the extrapolation velocity υe (τ). An event with a velocity error υ (τ) – υe (τ) is represented by an energy maximum either to the left or to the right of the Δτ = 0 line. The δτ(τ) trend can be picked and translated into a velocity trend as depicted in Figure 5.4-31b. This type of analysis has come to be called focusing analysis in the industry (Faye and Jeannaut, 1986). It has been used in some cases erroneously to estimate and update velocity-depth models used for depth migration. The method can only provide plausable velocity update within the framework of time migration.

Figure 5.4-32 is a CMP stack from offshore Texas. A 7000-ft portion (64 midpoints each with 48 offset traces) of the profile was used for migration velocity analysis. For computational efficiency, the data were windowed into 1024-ms time gates with 50 percent overlap. The image planes for one particular midpoint are shown in Figure 5.4-33. Different extrapolation velocities picked from a specified regional velocity function are used in each time gate. The velocity scan used in mapping is then carried out within a corridor around this function. Because different extrapolation velocities are used in successive segments, a given event appears at different values of τ in adjacent time segments.

The resulting velocity analysis for the central midpoint is shown in Figure 5.4-34. In conventional practice, to improve the quality of velocity picks, velocity analyses from a number of neighboring CMP gathers often are summed. Figure 5.4-34c shows the result of stacking velocity analysis for data from the six adjacent CMP gathers indicated in Figure 5.4-34a. For the migration-based method, the υ – τ planes corresponding to these gathers were summed. The result is shown in Figure 5.4-34b. The most obvious difference between the two results is the lack of shallow information in the migration-based υ – τ plane. This shortcoming is attributed to spatial aliasing and lack of long-offset data in the shallow time gate. The problem can be eliminated partly by increasing the length of the time gate used in the velocity analysis. With the shortcut time-windowing approach described above, the shallowest time segment did not include the large-offset data necessary for velocity resolution. Because the events have dip, the derived migration velocities are lower (by up to 4.5 percent) than the velocities derived from the stacking velocity analysis.

The velocity analysis described in this section does not handle lateral variations in velocity. It is based on a Fourier-transform domain formulation with only vertically varying velocity used in extrapolation. This method may be particularly efficient for the dip-corrected velocity estimate needed for time migration.

## Fowler's velocity-independent prestack migration

We now restate the underlying principle for migration velocity analysis:

Starting with the prestack volume of data P(y, h, t) in midpoint y, offset h and two-way event time t in the unmigrated position, create a velocity cube -- volume of data P(y, υ, τ) in midpoint y, migration velocity υ and two-way event time τ in the migrated position. Within the context of time migration, the output time variable τ is related to depth by way of vertical stretch:

$\tau ={\frac {2z}{\upsilon }}\$ Although the velocity cube can be created by means of some of the migration velocity analysis techniques described in this section, a variation of the method by Fowler (1984) is particularly efficient and elegant. First, we review Fowler's method to create the velocity cube. Refer to the moveout equation (5-1) and recall that stacking velocity υstk is dip dependent:

 $t^{2}=t_{0}^{2}+{\frac {4h^{2}}{\upsilon _{stk}^{2}}},$ (1)

where

 $\upsilon _{stk}={\frac {\upsilon _{DMO}}{\cos \phi }}.$ (2)

Use equation (5-11) to establish a relationship between the dip-dependent stacking velocities υstk and the dip-independent DMO velocities υDMO — velocities estimated from dip-moveout-corrected data:

 $\upsilon _{stk}={\frac {\upsilon _{DMO}}{\sqrt {1-\left({\frac {\upsilon _{DMO}k_{y}}{2\omega _{0}}}\right)^{2}}}}.$ (3)

This equation is the basis for Fowler mapping of υstk to υDMO. Note that the process is applied to data in the frequency-wavenumber domain. The Fowler mapping is then followed by constant-velocity Stolt mapping (Sections 4.1 and D.7) to map the DMO velocities υDMO to migration velocities υmig.

Stolt migration of the dip-moveout-corrected data volume in the Fourier transform domain P(ky0; υDMO) involves, first, mapping from ω0 to ωτ for a specific ky by using the dispersion relation of equation (4-24b) recast as

 $\omega _{0}={\sqrt {\omega _{\tau }^{2}+{\frac {\upsilon _{mig}^{2}k_{y}^{2}}{4}}}},$ (4)

where the relationship ωτ = υmigkz/2 is used. The output of this mapping P(kyτ; υmig) is then scaled by the quantity S given by equation (4-24c) recast as

 $S={\frac {\upsilon _{mig}}{2}}{\frac {\omega _{\tau }}{\sqrt {\omega _{\tau }^{2}+{\frac {\upsilon _{mig}^{2}k_{y}^{2}}{4}}}}}.$ (5)

Again, the relationship ωτ = υmigkz/2 is used to obtain equation 5 from equation (4-24c).

Figure 5.4-35 describes a flowchart for creating the migration velocity volume P(y, υmig, τ) by Fowler and Stolt mapping.

1. Start with data P(y, h, t) in coordinates of midpoint y, offset 2h and event time t in the unmigrated position, and create constant-velocity stack (CVS) volume P(y, υstk,tn) using a range of velocities υstk, where tn is the event time after constant-velocity normal moveout correction.
2. Apply 2-D Fourier transform to obtain the CVS cube P(ky, ustk) in the frequency-wavenumber domain, where ky and ω are the Fourier transform variables associated with the variables y and tn.
3. Sort the CVS volume P(ky, υstk) into a set of constant-velocity sections P(ky,ω; υstk).
4. Perform the Fowler mapping based on equation 3 on each of the velocity sections so as to obtain the DMO velocity volume P(ky0; υDMO).
5. Migrate each of the constant-velocity sections of the DMO velocity volume by performing the Stolt mapping based on equations 4 and 5 so as to obtain the migration velocity volume P(Ky, ωτ; υmig).
6. Apply 2-D inverse Fourier transform to obtain the migration velocity volume P(y, τ, υmig).

A variation of Fowler's sequence described above involves creating the CVS cube directly from DMO-corrected data.

1. Start with data P(y, h, t) in coordinates of midpoint y, offset 2h and event time t in the unmigrated position, and apply DMO correction followed by inverse NMO correction.
2. Create constant-velocity stack (CVS) volume P(y, υDMO,t0) using a range of velocities υDMO, where t0 is the event time after constant-velocity normal moveout correction (Figure 5.4-36a).
3. Sort the CVS volume P(y, υDMO,t0) into a set of cosntant-velocity stacked sections P(y,t0; υDMO).
4. Apply 2-D Fourier transform to obtain the CVS cube P(ky, υDMO0) in the frequency-wavenumber domain, where ky and ω are the Fourier transform variables associated with the variables y and t0.
5. Sort the CVS volume P(kyDMO0) into a set of constant-velocity sections P(ky0DMO).
6. Migrate each of the constant-velocity sections of the DMO velocity volume by performing the Stolt mapping based on equations 4 and 5 so as to obtain the migration velocity volume P(y, ωτ, υmig).
7. Apply 2-D inverse Fourier transform to obtain the migration velocity volume P(y, τ, υmig) (Figure 5.4-36b).

The migration velocity volume P(y, τ, υmig) shown in Figure 5.4-36b can be visualized and interpreted to derive a migration velocity field. For spatial consistency, velocity picking should be done on time slices from the migration velocity volume as shown in Figure 5.4-37. The resulting velocity strands shown in Figure 5.4-38a are interpolated to create the migration velocity field shown in Figure 5.4-38b. This velocity field then can be used to extract from the volume the section that corresponds to prestack time migration as shown in Figure 5.4-39. An enlarged view of this section is shown in Figure 5.4-40. Note the excellent imaging of the steeply dipping fault planes which conflict with the gently dipping strata.

Alternatively, the plane of (υmig) for each midpoint y can be inverse transformed to the plane of (h, τ) associated with the common-reflection-point gather derived from prestack time migration. This is then followed by conventional normal-moveout correction and stacking. The resulting section, again, represents the image from prestack time migration.

The data example shown in Figure 5.4-36 demonstrates the use of the migration velocity volume in deriving a high-fidelity image of fault planes from prestack time migration (Figure 5.4-40). The data example shown in Figure 5.4-41 demonstrates the use of the migration velocity volume in imaging steeply dipping event. The migration velocity volume was created using the procedure described above. Specifically, the DMO-corrected CMP gathers were first NMO-corrected using a range of constant velocities and a CVS volume was created. Next, each constant-velocity stacked section was migrated using the Stolt method and the constant velocity associated with that section. The resulting migration velocity volume is shown in Figure 5.4-41. Note the vertical variation in velocities on the end-on view of the volume that represents the plane of velocity versus time.

The migration velocity volume is interpreted by picking the primary velocity trend from selected time slices as shown in Figure 5.4-42a. Note the lateral variation in velocities which is captured by continuous picking along the midpoint axis. By interpolating the velocity strands resulting from the interpretation of selected time slices (Figure 5.4-42a), an rms velocity surface is generated. The picked velocity strands are shown in Figure 5.4-42b embedded within the migration velocity volume, and the rms velocity field is shown in Figure 5.4-43a as a color-coded surface extracted from within the migration velocity volume. A quality control of the rms velocity surface can be made by intersecting it with the cross-sections of the migration velocity volume at selected CMP locations as shown in Figure 5.4-43b.

The prestack time-migrated section is a by-product of the migration velocity analysis described here. Specifically, the image surface associated with the prestack time migration is obtained by sculpting the amplitudes from within the migration velocity volume over the rms velocity surface as shown in Figure 5.4-44a. The conventional 2-D display of the prestack time-migrated section is then created by simply collapsing the sculpted image surface onto a 2-D plane (Figure 5.4-44b). A close-up view of the prestack time-migrated section shows accurate imaging of the steeply dipping event (Figure 5.4-45).

## Exercises

Exercise 1
Consider the application of DMO correction to data referred to a floating datum represented by a smooth form of an irregular topographic surface and data referenced to a flat datum below. Which DMO correction would have more effect on the data?
Exercise 2
Refer to Fowler's velocity-independent prestack migration technique for migration velocity analysis. Suppose you have transformed the prestack data from offset to velocity space using 30 constant velocity values from 2000 m/s to 4900 m/s using an increment of 100 m/s. Can you create additional constant-velocity panels such that the increment is 50 m/s by poststack migration rather than prestack migration or trace interpolation? If so, what velocity would you use for poststack migration to create the panel for 3050-m/s velocity.
Exercise 3
Derive Bancroft's equivalent offset equation (E-72b) from the nonzero-offset traveltime equation (E-67).
Exercise 4
Suppose you have two events with conflicting dips of the same magnitude, but in opposite directions, associated with a reflector within a fault block and the fault plane itself. Can these two events be distinguished on a velocity spectrum computed from a CMP gather before DMO correction at a location just above the surface point where the two events intersect one another?

## Migration velocity analysis

The idea of a velocity analysis that is based on differential solutions of the scalar wave equation first was introduced by Doherty and Claerbout . They used the 15-degree finite-difference migration algorithm and worked with single CMP gathers. Gonzalez-Serrano and Claerbout  later extended the wave equation velocity analysis to slant-midpoint coordinates and worked with linearly moveout-corrected CMP gathers. The method discussed here  operates in the Fourier transform domain using the exact form of the double-square-root (DSR) operator. Mathematical details of this method are found in Section E.7. Figure 5.4-23 summarizes the main computational steps involved in this migration velocity analysis based on wavefield extrapolation.

1. Starting with the prestack data in midpoint y, offset h, and two-way event time t in the unmigrated position, represented by the wavefield P(y, h, τ = 0, t) at the surface τ = 0, perform 3-D Fourier transform. The variable τ is associated with the direction of wave extrapolation, and is related to depth z by τ = 2z/v, where v is the medium velocity.
2. Specify an extrapolation velocity function that only varies vertically, v(τ) and apply the extrapolation operator exp(−iωDSRτ/2) to compute the extrapolated wavefield in the transform domain P(ky, kh, τ, ω) from the surface wavefield in the transform domain P(ky, kh, τ = 0, ω).
3. To obtain the zero-offset image, sum over the offset wavenumber, and thus obtain P(ky, h = 0, τ, ω).
4. Apply 2-D inverse Fourier transform to obtain the zero-offset image P(y, h = 0, τ, t). The image below a midpoint y is contained in the t − τ plane.
5. Perform mapping of the variables as described in Section E.6 from τ to v. The velocity information is given by the envelope of the velocity volume of data P(y, h = 0, τ = t, v).

We now demonstrate the procedure outlined in Figure 5.4-23 using a synthetic data set. Figure 5.4-24 shows two common-offset sections over a number of point scatterers buried in a constant-velocity earth, where v = 3000 m/s. Using a constant velocity for extrapolation, ve = 3000 m/s, the t − τ image plane was produced for each midpoint. Two such planes corresponding to midpoints 1 and 5 denoted in Figure 5.4-24 are shown in Figure 5.4-25. The v − τ planes (Figure 5.4-26) then were generated from the t − τ image planes by the mapping procedure described in Section E.7. Peak amplitudes for all events occur at the correct medium velocity (3000 m/s). We expect the diffraction events in Figure 5.4-23 to migrate to the apexes beneath midpoint 1, where the point scatterers are located. Note that in Figure 5.4-25, almost all the energy is in the image plane corresponding to midpoint 1; just five midpoints away, at midpoint 5, the migrated energy is very low.

How do we interpret the t − τ image planes? If we used the true medium velocity in downward extrapolation, then, according to the imaging principle, we would see all the events along the diagonal τ = t, the image line, on the image plane. This happens in Figure 5.4-25, because a 3000-m/s extrapolation velocity was used, which is just the velocity used in generating the model in Figure 5.4-23. Any displacement of peak energy from the t = τ image line means that the velocity value used for downward extrapolation differs from that of the event. This displacement is also the basis for mapping from the t − τ image plane to the v − τ plane by equation (E-77).

 $Y={\frac {vk_{y}}{2\omega }}$ (E-77a)

 $H={\frac {vk_{h}}{2\omega }}.$ (E-77b)

This mapping is investigated further with the modeled data set shown in Figure 5.4-27, in which velocity increases with depth. In Figure 5.4-28b, note that the top and middle events fall to the left of the image line suggesting that the velocity used in extrapolation (ve = 3000 m/s) is greater than the velocities associated with these events. The bottom event falls on the image line, implying that its velocity is nearly the same as the extrapolation velocity. These observations are confirmed in the corresponding v − τ planes in Figure 5.4-29. While true stacking velocity values for the three events are 2700, 2850, and 3000 m/s, the velocities interpreted from Figure 5.4-29b are about 2500, 2800, and 3000 m/s. Thus, the migration-based velocity estimate for the shallow event is in error by approximately 8 percent.

To determine the reason for the velocity error, we will consider a migration-based velocity analysis of our synthetic data example that does not involve the approximate mapping step. Figure 5.3-30a shows a CMP gather from midpoint 1 in the zero-dip region of the depth-variable velocity model associated with the constant-offset sections in Figure 5.4-27. The migration velocity analysis on this gather (Figure 5.4-30b) was done by extrapolating the surface wavefield P(kh, ω, τ = 0) repeatedly with different constant velocities in steps of Δτ = Δt (the sampling rate). The zero-offset trace from each attempt was collected after this effort, abandoning the rest of the migrated CMP gather.

Interpretation of the velocity analysis in Figure 5.4-30b reveals correct stacking velocities for the three events, including the shallowest. Clearly the error observed in Figure 5.4-29 is attributable to the mapping (E-100). Note that the error does not occur because of depth variability of velocity, but instead, because the single extrapolation velocity used differed from the medium velocity. The conventional velocity analysis for midpoint 1 of this model data set is shown in Figure 5.4-30c for comparison. Note the familiar NMO stretching that is apparent in the shallow event. In other respects, both the results (Figures 5.4-30b and 5.4-30c) are comparable.

 $\tau \,v_{e}^{2}=t\,v^{2}.$ (E-100) Figure 5.4-33  Image planes associated with the center midpoint from the zone of interest in Figure 5.4-32.

The departure of an event on the t − τ image plane from the t = τ image line is measured by the quantity Δτ as depicted in Figure 5.4-31a. In some practical implementations, the t − τ image plane is mapped onto the plane of Δτ versus τ as depicted in Figure 5.4-31c to determine the rms velocity v(τ) for time migration from the extrapolation velocity ve(τ). An event with a velocity error v(τ) − ve(τ) is represented by an energy maximum either to the left or to the right of the Δτ = 0 line. The δτ(τ) trend can be picked and translated into a velocity trend as depicted in Figure 5.4-31b. This type of analysis has come to be called focusing analysis in the industry . It has been used in some cases erroneously to estimate and update velocity-depth models used for depth migration. The method can only provide plausable velocity update within the framework of time migration.

Figure 5.4-32 is a CMP stack from offshore Texas. A 7000-ft portion (64 midpoints each with 48 offset traces) of the profile was used for migration velocity analysis. For computational efficiency, the data were windowed into 1024-ms time gates with 50 percent overlap. The image planes for one particular midpoint are shown in Figure 5.4-33. Different extrapolation velocities picked from a specified regional velocity function are used in each time gate. The velocity scan used in mapping is then carried out within a corridor around this function. Because different extrapolation velocities are used in successive segments, a given event appears at different values of τ in adjacent time segments.

The resulting velocity analysis for the central midpoint is shown in Figure 5.4-34. In conventional practice, to improve the quality of velocity picks, velocity analyses from a number of neighboring CMP gathers often are summed. Figure 5.4-34c shows the result of stacking velocity analysis for data from the six adjacent CMP gathers indicated in Figure 5.4-34a. For the migration-based method, the v − τ planes corresponding to these gathers were summed. The result is shown in Figure 5.4-34b. The most obvious difference between the two results is the lack of shallow information in the migration-based v − τ plane. This shortcoming is attributed to spatial aliasing and lack of long-offset data in the shallow time gate. The problem can be eliminated partly by increasing the length of the time gate used in the velocity analysis. With the shortcut time-windowing approach described above, the shallowest time segment did not include the large-offset data necessary for velocity resolution. Because the events have dip, the derived migration velocities are lower (by up to 4.5 percent) than the velocities derived from the stacking velocity analysis. Figure 5.4-34  (a) A portion of the CMP stacked section in Figure 5.4-32, (b) the velocity spectrum based on the procedure in Figure 5.4-23, followed by averaging over six midpoints; (c) the conventional velocity analysis as discussed in velocity analysis followed by averaging over six midpoints.

The velocity analysis described in this section does not handle lateral variations in velocity. It is based on a Fourier-transform domain formulation with only vertically varying velocity used in extrapolation. This method may be particularly efficient for the dip-corrected velocity estimate needed for time migration.