I 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,
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).
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
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 X2. 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.
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
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
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
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
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.
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.
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.
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 15 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.
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).
Figure 17 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 19 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 20 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.
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.