biondo@sep.stanford.edu

## ABSTRACTI derive an unconditionally stable implicit finite-difference operator that corrects the constant-velocity phase shift operator for lateral velocity variations. My method is based on the Fourier-Finite Difference (FFD) method first proposed by Ristow and Ruhl 1994. Contrary to previous results, my correction operator is stable even when the reference velocity is higher than the medium velocity. Because of this additional capability, after the correction step I can apply a frequency-dependent interpolation that significantly reduces: the residual phase error after correction, the frequency dispersion caused by the discretization of the Laplacian operator, and the azimuthal anisotropy caused by splitting. Tests on zero-offset data from the SEG-EAGE salt data set show that the proposed method improves the imaging of a fault reflection with respect to a similar interpolation scheme that employs extended split-step to adapt to lateral velocity variations. |

As 3-D prestack wave-equation imaging is becoming practically possible Biondi and Palacharla (1996); Mosher et al. (1997), we need robust, efficient, and accurate methods to downward continue 3-D wavefields. In particular, wide-angle methods are crucial for prestack imaging because at least one of the paths connecting the image point in the subsurface to the source/receiver locations at the surface is likely to propagate at a wide angle.

Fourier methods, such as phase shift Gazdag (1978), handle wide-angle propagation efficiently and accurately, but only for vertically layered media. In contrast, finite-difference methods can easily handle lateral velocity variations, but are not efficient for wide-angle propagation. A natural strategy is thus to combine a Fourier method with a finite-difference method to derive an extrapolation method that enjoys the strengths of both. This is not a new idea, and indeed both methods that were first proposed to adapt Fourier methods to lateral velocity variations, Phase Shift Plus Interpolation (PSPI) Gazdag and Sguazzero (1984) and split-step Stoffa et al. (1990), can be interpreted as being zero-order finite-difference corrections to a phase shift extrapolator. Ristow and Ruhl 1994 first proposed a genuinely finite-difference correction to phase shift, which they dubbed Fourier-Finite Difference (FFD), that employs implicit finite differences Claerbout (1985) to handle lateral velocity variations. Pseudo-screen propagators Jin et al. (1998), generalized-screen propagators Le Rousseau and de Hoop (1998), and local Born-Fourier migration Huang et al. (1999) are related methods that have been recently proposed in the literature.

In the first part of this paper, I show that the FFD correction is more accurate than the implicit finite-difference implementation of the first order Born-Fourier and pseudo-screen corrections. Because the computational complexity of the three methods is comparable, the FFD correction is more attractive than the others, and it is the focus of my paper. Unfortunately, my first attempts at applying the original FFD method failed because of numerical instability arisen when the velocity model had sharp discontinuities, such as the one encountered at the boundaries of a salt body. Stability is a necessary condition for a migration method to be practically useful.

The stable FFD correction that
I present in this paper
overcomes the instability problems related
to the original FFD method.
To derive a stable version of the FFD correction
I adapted the bullet-proofing theory developed
by Godfrey et al. 1979
and Brown 1979
for the 45-degree equation.
The bullet-proofed FFD correction
is unconditionally stable for arbitrary variations
in the medium velocity *and*
in the reference velocity.
Further, I demonstrate that it is unconditionally stable
when the medium velocity is either higher *or*
lower than the reference velocity.
This is a useful result,
and differs with
a statement in Ristow and Ruhl's paper,
that asserts their method to be unstable when the
medium velocity is lower than the reference velocity.

The stability of the new FFD correction,
even when the reference velocity is higher than the medium velocity
and it is laterally varying,
makes it a suitable building block for the construction
of a new wide-angle downward continuation
algorithm that is efficient and accurate in 3-D.
At each depth step, the wavefield is propagated
with *n* reference velocities using phase shift,
where *n* is determined according to the
maximum propagation angle that we need to propagate accurately.
Then the *n* reference wavefields are combined to create two wavefields:
one for which the reference velocity is lower than the medium velocity,
the other one for which the reference velocity is higher than the medium
velocity.
A stable FFD correction is applied to both wavefields,
and the corrected wavefields are linearly interpolated
with frequency-dependent weights.
The frequency-dependent interpolation enables a
significant reduction of the frequency dispersion introduced
by the discretization of the Laplacian operator
in the implicit finite difference step.

In 3-D, the FFD corrections can be efficiently applied by splitting, or possibly by helix-transform methods Rickett (2000). However, the proposed algorithm suffers much less from azimuthal anisotropy caused by splitting than the original FFD method. The phase variations as a function of azimuth have opposite behavior when the differences between the reference velocity and medium velocity have opposite signs. Therefore, these phase variations tend to cancel each other when the two wavefields are interpolated after the correction.

Because Fourier-finite differences methods and interpolation are both fundamental components of the new method, I will refer to it as Fourier-Finite Difference Plus Interpolation method, or FFDPI.

Pseudo-screen vs. FFD correction
The simplest way to derive the
pseudo-screen, or local Born-Fourier, correction
is to expand
the square root of the one-way wave equation
in a Taylor series and stop at the first order.
The wavefield downward continued at depth is then computed from the wavefield at depth *z*
using the following approximation,

(1) |

(2) |

(3) |

(4) |

(5) |

Figure 1 demonstrates the accuracy improvement
gained by including the second term in equation (5).
It compares the phase curves
obtained after applying the first term in
equation (5) (split step)
and after applying both terms (pseudo screen).
The medium velocity *v* is equal to 2 km/s,
and two reference velocities are assumed:
one 10% lower than the medium velocity (1.8 km/s),
the other one 10% higher than the medium velocity (2.2 km/s).

Figure 1

FFD correction
The FFD correction achieves better accuracy than
the pseudo-screen correction
because it is based on a direct expansion of
the difference between the
square root evaluated at the medium velocity *v*
and the square root evaluated at the reference velocity ,instead of being based on the expansion of the square root
around the reference velocity.
The downward continued wavefield is approximated as

(6) |

(7) |

(8) |

Figure 2

Figures 3-5
show the impulse responses associated with the phase curves
shown in Figure 2.
The maximum frequency in the data is 42 Hz and the
spatial sampling is 10 m in both directions.
Figure 3 shows the exact impulse
response for the medium velocity equal to 2 km/s.
Figure 4 shows the impulse response
with reference velocity equal to 1.8 km/s and FFD correction.
Figure 5 shows the impulse response
with reference velocity equal to 2.2 km/s and FFD correction.
Notice the frequency dispersion in the image obtained
applying the FFD correction.
These artifacts are caused by the discretization errors
of the horizontal second derivatives in *X ^{2}*.
The phase curves shown in Figure 2
neglect this approximation,
and thus they represent the effective phase shift for zero-frequency data.
Also notice that the frequency dispersion
is in the opposite directions for opposite signs of
the velocity correction.

Figure 3

Figure 4

Figure 5

Stable FFD correction An implicit finite-difference implementation using a Crank-Nicolson scheme of the FFD correction as expressed in equation (8) is stable for smooth velocity variations. But numerical instability may develop when there are sharp discontinuities in the velocity field. An example of this situation is shown in Figures 6 and 7. The slowness function (Figure 6) has a sharp negative step, and a random behavior within the low slowness region. The impulse response computed applying the original FFD correction is shown in Figure 7, and clearly illustrates the problem. Notice that the image was clipped at the 70th percentile before plotting, to allow the ``shadow'' of the familiar circular impulse response to be visible in the plot.

In constant velocity the correction operator is unitary (all-pass filter) because its eigenvalues have zero imaginary part. Numerical instability originates when variations in the velocities terms multiplying the second derivative [ and ]causes the imaginary part to become different from zero. To assure that this does not happen we first rewrite equation (8) as

(9) |

(10) |

(11) |

The matrix
is now guaranteed to have
real eigenvalues.
Because both *I*
and
are normal matrices that can be diagonalized by the
same similarity transformation
Brown (1979),
also the matrix
is now guaranteed to have
real eigenvalues.

In matrix notation equation (9) can be rewritten as

(12) |

(13) |

(14) |

(15) |

(16) |

The reference velocity and the medium velocity *v*
can be interchanged at will in the previous development
without changing the stability conditions.
Therefore, the stable FFD correction is not only
stable in presence of sharp discontinuities
in the medium velocity, but also
in presence of sharp discontinuities
in the reference velocity.
This property is exploited in the next section
to design an efficient and accurate interpolation scheme.

Equation (14) can be solved using a Crank-Nicolson scheme and the wavefield at depth computed as

(17) |

Figure 8, shows the same impulse response as in Figure 7, but computed using the stable FFD correction. In this case no numerical instability is encountered and the wavefield propagates without problems through the region with random slowness perturbations.

In 3-D the stable FFD correction could be applied using a 2-D version of the Crank-Nicolson method summarized in equation (17). However, the solution of the linear system implied by equation (17) would be expensive. A splitting algorithm Jakubowicz and Levin (1983) can be employed to reduce computational costs. The stability of a splitting algorithm derives directly from the analysis above, since it consists of the successive application of the FFD correction along the two horizontal coordinate axes. In the next section I discuss how the capability of using both positive and negative velocity corrections yields a significant improvement in accuracy of the splitting algorithm.

Boundary conditions
A necessary component of deriving a stable
downward continuation scheme
is to define stable boundary conditions.
It is also desirable for the boundaries to be absorbing.
Following Clayton and Engquist 1977,
and Rothman and Thorson
1982,
this goal can be easily accomplished by
changing the values at the edges of the diagonal
of *T* in equation (10)
from 2 to ,where
;that is substituting *T*
with

(18) |

Figure 6

Figure 7

Figure 8

The Fourier-Finite Difference Plus Interpolation algorithm The stable FFD correction developed in the previous section has the right characteristics to be used as the main building block of an efficient and accurate wide-angle continuation algorithm. To achieve accuracy we can interpolate between wavefields that have been phase shifted using several reference velocities and corrected using the stable FFD method. In theory, arbitrary accuracy can be achieved by increasing the number of reference velocities. The structure of the algorithm is similar to the PSPI method, except that a wide-angle correction (FFD) is employed instead of a narrow-angle one (vertical shift). This change reduces the errors over the whole range of propagation angles.

Two results reached in the previous section
are important for the definition of a stable and accurate
interpolation scheme.
First, the application of the FFD correction is stable independently
from the sign of the velocity correction to be applied,
as long as it is constant within the same correction step.
This result enables a *linear interpolation*
between a wavefield corresponding to reference velocities
lower than the medium velocity
and a wavefield corresponding to reference velocities
higher than the medium velocity.
Previously, because of the requirement for the
reference velocity to be lower than the medium velocity,
only a *nearest-neighborhood interpolation* was possible.
Second, the reference velocity can vary at will laterally,
without compromising the stability of the method.
Because of these results,
computations can be saved
by applying the FFD correction
only twice at each depth step.
Once starting from a wavefield constructed from all the reference wavefields
computed using a reference velocity lower than the medium velocity.
The second time starting from
a wavefield constructed from all the reference wavefields
computed using a reference velocity higher than the medium velocity.

The algorithm outlined above can be described in more details as the sequence of the following steps.

- 1.
- Downward continue the data using phase shift for reference velocities ,and compute the reference wavefields as
(19) - 2.
- Define two reference velocity functions and ,that at every point are respectively equal to
the reference velocity that is just lower and just higher
than the medium velocity;
that is,
(20) (21) - 3.
- Extract from the reference wavefields two
wavefields corresponding to and and correct them using the stable FFD method:
(22) (23) - 4.
- Linearly interpolate the two corrected wavefield as
(24) (25) *k*_{m}. When the second derivatives are computed using the first order approximation*T*_{b}in equation (18), is given as a function of*k*_{m}by the following expression,(26)

It is important to notice that the stability analysis developed in the previous section strictly applies only to the simple FFD correction, not to its combination with an interpolation scheme like FFDPI. In theory, instability can still develop when using FFDPI, as it does for PSPI Margrave and Ferguson (1999). However, in my tests I have not encountered it yet!

FFDPI error analysis
Figure 9 shows the effect of interpolation on the
phase curves after interpolation.
The phase function was interpolated after
the application of the FFD correction starting from a reference
velocity lower than the medium velocity
and a reference velocity higher than the medium velocity.
As in Figure 2, the medium velocity *v* is equal to 2 km/s,
and two reference velocities are assumed:
one 10% lower than the medium velocity (1.8 km/s),
the other one 10% higher than the medium velocity (2.2 km/s).
The interpolation weights were computed using
equation (26) at an angle of 60 degrees.

Figure 10 shows the phase error measured
at 50 degrees as a function of the medium velocity *v*
as it spans the range between
the lower reference velocity (1.8 km/s)
and the higher reference velocity (2.2 km/s).
It compares the phase errors at three different frequencies
(0 Hz, 30 Hz and 60 Hz)
when the interpolation weights are computed assuming
an exact second derivative operator
(i.e. weights are constant with frequency)
and the first order approximation of the second derivative operator
(i.e. weights vary with frequency).
The absolute value of the maximum error suffered
when frequency-dependent weights are used
(.1% at 60 Hz and 1.95 km/s)
is very small,
and almost a factor of 8 smaller than when
frequency-independent weights are used.

Figure 9

Figure 10

Figures 11 and 12 show the impulse responses corresponding to the phase plots shown in Figures 9 and 10 and should be compared with the impulse responses shown in Figures 3-5. The impulse response shown in Figure 11 was computed using frequency-independent interpolation weights. While it is much closer to the exact impulse response (Figure 3) than either of the impulse responses obtained using a simple FFD correction (Figures 4 and 5), it shows some frequency dispersion. The higher frequencies are imaged inside the semicircle. The frequency dispersion is greatly reduced when the frequency-dependent interpolation weights are used, as demonstrated in Figure 12, and predicted by the plots in Figure 10.

A wavefield interpolation scheme
similar to the one described above
could be used in conjunction with the split-step correction
instead of the FFD correction.
With a fixed number of reference velocities,
using split step instead of FFD would reduce the computational cost.
It is thus useful to compare the phase errors of the
two competing methods,
to analyze the improvement in accuracy gained using the more accurate
but more expensive correction method.
Figure 13 compares the phase errors measured
at 50 degrees as a function of the medium velocity *v*
as it spans the range between
the lower reference velocity (1.8 km/s)
and the higher reference velocity (2.2 km/s).
The maximum error suffered by the method based on the
split-step correction is about 8 times larger than the error suffered
when the FFD correction is applied.

Figure 11

Figure 12

Figure 13

Azimuthal anisotropy A recurring problem that hampers the application of implicit finite-difference methods to 3-D wave extrapolation is the azimuthal anisotropy associated with splitting Jakubowicz and Levin (1983). Of course, this problem affects also the FFD correction applied by splitting 1998. One potentially attractive way of solving this problem is using helical boundary conditions, as Rickett discusses in an article in this report 2000. However, I will show that for the FFDPI algorithm the azimuthal anisotropy is greatly reduced without recurring to sophisticated linear solvers.

Figure 14 compares relative phase errors
as a function of the azimuth measured for a propagation
angle of 53 degrees.
As in the previous figures,
the medium velocity *v* is equal to 2 km/s,
and two reference velocities are assumed:
one 10% lower than the medium velocity (1.8 km/s),
the other one 10% higher than the medium velocity (2.2 km/s).
The frequency-dependent interpolation weights were computed
to zero the phase error along the inline/crossline
directions for the same propagation angle.
The plots show the phase errors at two frequencies (0 Hz and 60 Hz)
for the FFDPI algorithm,
the FFD correction starting from
the lower reference velocity,
and
the FFD correction starting from
the higher reference velocity.
Notice that for both the simple FFD correction cases
the azimuthal anisotropy decreases as the frequency increases,
though the average phase error increase as well.
But the crucial, and useful,
feature of the phase errors function
for the FFD corrections is that the azimuthal variations
are in the opposite direction
when the differences between the reference velocity
and medium velocity have opposite sign.
Therefore, the phase error of the interpolation method
is much lower than the error of either of
the simple FFD corrections.
At higher frequencies (60 Hz) the impulse response
of FFDPI is almost perfectly isotropic.

The theoretical analysis is confirmed by the characteristics of the impulse responses. Figure 15 shows the depth slice of three impulse responses superimposed to each other. The outermost circular event corresponds to the FFD correction starting from a reference velocity of 2.2 km/s. The middle event corresponds to the exact impulse response with the medium velocity of 2 km/s. The innermost event corresponds to the FFD corrections starting from a reference velocity of 1.8 km/s. As predicted by the plots shown in Figure 14, the azimuthal anisotropy is frequency dependent and the frequency dispersion is smaller for azimuths oriented at 45 degrees with respect to the coordinate axes.

Figure 16 is the merge of two impulse responses along the inline direction. For negative values of the in-line coordinate the plot shows the depths slice for the exact impulse response. For positive values of the in-line coordinate the plot shows the depths slice for the impulse response obtained by FFDPI. It is evident that the result of the interpolation scheme is much less affected by azimuthal anisotropy and frequency dispersion than the results of the two simple FFD correction showed in Figure 15.

Figure 14

Aniso-rist
Depth slices through impulse responses:
1) innermost event corresponds to the FFD corrections starting
from a reference velocity of 1.8 km/s,
2) middle event corresponds to the exact impulse response
with the medium velocity of 2 km/s,
3) outermost event corresponds to the FFD corrections starting
from a reference velocity of 1.8 km/s.
Figure 15 |

Figure 16

Zero-offset migration of the SEG-EAGE salt data set To test both the stability and the accuracy of the FFDPI algorithm, I migrated the zero-offset data from the SEG/EAGE salt data set Aminzadeh et al. (1996). The data set is a good test for the stability of the FFDPI algorithm because the velocity model has sharp discontinuities caused by the salt body. Furthermore, because of a low-velocity region intended to model subsalt overpressure, several depth slices have a wide range of velocities. Figure 17 shows one of these depth slices. In the plot the salt velocity is clipped, thus the scale-bar on the side represents the range of velocities within the sediments. There is almost a factor of two between the slow velocity sediments in the `overpressure zone' in the middle, and the faster sediments at the edges.

The migration algorithm does not need to handle accurately lateral velocity variations to image the reflectors above the salt. While the reflectors below the salt cannot be imaged by simple zero-offset migration. A deep fault (between depths of 2 km and 3 km) and away from the salt body is one of the few reflectors that is well suited to test the accuracy of a zero-offset migration. In previous reports I noted that six reference velocities are needed to image this particular fault by common-azimuth migration employing a split-step correction Biondi (1999a,b). Figure 18 shows an in-line section of the migrated cube that cuts across the fault of interest. It was obtained using the FFDPI algorithm. Four reference velocities were used at each depth step. Figures 19 and 20 compare the zooms around the fault of interest. The light (yellow) line superimposed onto both plots represents the correct fault position, as picked from the velocity model. Figure 19 shows the results obtained by using the interpolation algorithm described above, but in conjunction with the split-step correction instead of the stable FFD correction. The fault is undermigrated and thus it is imaged too shallow, and the sediment terminations are not well focused. The fault is better positioned and the image is better focused when the FFDPI algorithm is used to migrate the data (Figure 20).

Vel_spk_z2100
Depth slice of the velocity model at z=2.1 km.
The salt velocity was clipped, thus the scale bar on the side
shows the range of velocities in the sediments.
Figure 17 |

Figure 18

salt-spl
Window of the same in-line section shown in
Figure 18, but obtained by
use of split step to correct for lateral velocity variations.
The light (yellow) line superimposed onto the plot represents the correct fault
position, as picked from the velocity model.
Notice the misplacement of the fault.
Figure 19 |

salt-rist
Window of the same in-line section shown in
Figure 18 and obtained by
use of FFDPI.
The light (yellow) line superimposed onto the plot represents the correct fault
position, as picked from the velocity model.
Notice the better placement of the fault.
Figure 20 |

Conclusions The combination of Fourier methods' accuracy for wide-angle propagation with implicit finite differences' flexibility for modeling lateral velocity variations yields accurate and efficient downward-propagation methods. The FFD corrections is the most attractive among several methods proposed in the literature to correct constant-velocity phase shift for lateral velocity variations. However, the correction operator originally presented by Ristow and Ruhl 1994 can be unstable in presence of sharp discontinuities in the velocity function. In this paper I present and successfully test a stable version of the FFD correction.

Using the stable FFD correction as a building block, I derive an accurate and stable wide-angle migration (Fourier-Finite Difference Plus Interpolation). The FFDPI algorithm is based on the interpolation of two wavefields corrected using the FFD method with opposite signs of the velocity perturbations. This interpolation step compensates for both the azimuthal anisotropy and the frequency-dispersion of the simple FFD corrections. Therefore the FFDPI algorithm achieves high accuracy, as demonstrated by the migration example of the SEG-EAGE salt data set. The accuracy and the cost of FFDPI algorithm can be easily controlled by setting the number of reference velocities.

4/27/2000