Abstract
Three-dimensional wave-equation migration techniques are still quite expensive because of the huge matrices that need to be inverted. Several techniques have been proposed to reduce this cost by splitting the full 3D problem into a sequence of 2D problems. We compare the performance of splitting techniques for stable 3D Fourier finite-difference (FFD) migration techniques in terms of image quality and computational cost. The FFD methods are complex Padé FFD and FFD plus interpolation, and the compared splitting techniques are two- and four-way splitting as well as alternating four-way splitting, that is, splitting into the coordinate directions at one depth and the diagonal directions at the next depth level. From numerical examples in homogeneous and inhomogeneous media, we conclude that, though theoretically less accurate, alternate four-way splitting yields results of comparable quality as full four-way splitting at the cost of two-way splitting.
1. Introduction
Because of its superiority in areas of complex geology, wave-equation migration is substituting Kirchhoff migration in practice. However, while Kirchhoff migration counts on more than 30 years of technological development, wave-equation migration methods still need to be improved in various aspects. One of these aspects is the efficient implementation of three-dimensional wave-equation migration.
The application of a three-dimensional wave-equation migration technique adds the problem of computational cost to those of stability and precision of the chosen migration algorithm. To speed up migration techniques like finite-difference (FD) [1] or Fourier finite-difference (FFD) migration [2], a technique known as splitting is frequently used. In this context, splitting means the separation of a single-step 3D migration into two 2D passes within planes parallel to the horizontal coordinate axes, usually the inline and crossline directions [3].
When the splitting is applied to the implicit FD migration operator in such a way that the resulting equations are solved alternatingly in the inline and crossline directions, the resulting FD scheme is known as an alternating-direction-implicit (ADI) scheme. This procedure has the drawback of being incorrect for strongly dipping reflectors, resulting in large positioning errors for this type of reflectors when the dip direction is away from the coordinate directions and thus outside the migration planes. This imprecision leads to numerical anisotropy, that is, a migration operator that acts quite differently in different directions.
To improve this behaviour while retaining the advantages of a rather low computation cost, different procedures have been proposed over the years. Ristow ([4], see also [5]) proposed to perform, in addition to the 2D migration in the coordinate planes, also 2D migrations in the diagonal directions between the coordinate axes. Kitchenside [6] used phase-shift migration plus an additional FD propagation step of the residual field to reduce the splitting error. Graves and Clayton [7] proposed the implementation of a phase-correction operator using FD and incorporating a damping function to guarantee the stability of the 3D FD migration scheme.
Inverting the idea of Kitchenside [6], who propagated the field using phase shift and the residual using FD, Li [8] proposed to use conventional FD migration plus a residual field correction by phase shift to improve the migrated image quality. Without any need to modify the conventional 3D FD migration, the Li correction adds a phase-shift filter at certain steps of the downward extrapolation. This technique corrects not only for the splitting error but also for the positioning error of steeply dipping reflectors.
Collino and Joly [9] solved a family of new 3D one-way wave equations by the ADI method. These equations significantly reduce the numerical anisotropy, but are approximately four times as expensive as conventional two-way splitting. Wang [10] developed an alternative method to improve the precision of the FD solution of the one-way wave equation. To guarantee stability and efficiency, he keeps the implicit FD scheme and the alternation of directions, but interpolates between the ADI solution and the wavefield before each step of extrapolation. He calls the resulting method ADI plus interpolation (ADIPI). As a drawback, this ADIPI can produce instabilities in the presence of strong lateral velocity variations.
Zhou and McMechan [11] proposed a one-way wave equation that can be expressed as a system of differential equations of first and second order [12] and factored into a product of two one-dimensional terms corresponding to the lateral directions [7, 13]. A big asset of the method is that the conventional FD extrapolation can be used with very little modification. In this way, the efficiency of conventional splitting is preserved, without adding the necessity of any error compensation. However, the method is also instable for strong lateral velocity contrasts and needs rather heavy model smoothing.
Biondi [14] showed that FFD migration is more precise than other methods that use implicit finite differences like pseudoscreen propagators [15] and high-angle screen propagators [16]. Given that the computational complexity of all three methods is approximately the same, FFD migration is more attractive than the others. Unfortunately, when conventional FFD migration is applied in the presence of strong velocity contrasts, it can generate numerical instabilities, too.
To overcome the problem of instabilities in models with strong lateral velocity contrasts, Biondi [14] presented a correction to the FFD method that avoids stability problems. To derive it, he adapted a theory of Godfrey et al. [17] and Brown [18], which improves the stability of the equation. The corrected FFD method is unconditionally stable for arbitrary velocity variations, as much in the velocity model as in the reference velocity. Particularly, and differently from conventional FFD migration, it is unconditionally stable even if the reference velocity is smaller than the model velocity. This new property allows for the application of the interpolation technique, conventionally used to improve phase-shift and split-step migration [19] but impossible in FFD migration, because it needs propagation with a larger and a smaller reference velocity. The resulting migration technique is called FFD plus interpolation, or shortly FFDPI.
Another computationally less expensive method to stabilize FFD migration in the presence of strong lateral velocity contrasts was proposed by Amazonas et al. [20]. It substitutes the real Padé approximation [21] used in the derivation of FFD migration [2] by its complex version [22]. In this way, the incorrect treatment of near horizontal and slightly evanescent waves of the real Páde approximation is improved, leading to a more stable FFD algorithm, shortly referred to as complex Padé FFD (CPFFD) migration.
In this work, we study possibilities of efficiently implementing these stable FFD migration techniques in 3D. We implemented and compared splitting techniques for FFDPI [14] and CPFFD [20] migration. Our numerical tests indicate that a very robust, highly efficient, and satisfactorily accurate method is alternate four-way splitting, that is, splitting into the coordinate directions at one extrapolation step and into the diagonal directions at the next step.
2. Theoretical Background
2.1. The One-Way Wave Equation
The one-way wave equation [23] can be derived starting from the scalar wave equation, which for a homogeneous medium is given by where is the scalar wavefield and is the spatially varying wave velocity. For moderately varying media, where velocity derivatives can be neglected, Fourier transform in time and horizontal coordinates and allows to represent (1) as Equation (2) can be factorized into The two differential operators in (3) represent, when taken alone, one-way wave equations that describe up- and downgoing waves. For migration, the one-way wave equation of interest is the one describing downgoing waves, that is, Inverse Fourier transform in the horizontal wavenumbers and yields then formally The actual restrictions that apply to (5) in inhomogeneous media are much less severe than the above derivation indicates. Of course, for the formal representation (5) to make practical sense, the square root of the differential operator needs to be approximated in terms of numerically executable operations.
2.2. Expansion of the Square Root
A well-used possibility for the approximation of the square root in the one-way wave equation (5) in terms of numerically executable operations is an expansion into a Padé series [21]: where the Padé coefficients are
This approximation is used in most practical FD migration schemes. Depending on the number of terms used in the expansion, it gives rise to the so-called , , or migrations.
However, when the interest is on accurate imaging up to very high propagation angles, approximation (6) has a drawback, because its validity is limited to . For , the approximation breaks down abruptly, because the left side of (6) is imaginary, while its right side remains real. Thus, for propagation angles close to , where the argument of (6) becomes close to −1, (6) becomes invalid, which causes instabilities when using the real approximation (6) for migration in models with strong lateral velocity contrasts.
To overcome this problem, Milinazzo et al. [22] proposed to rotate the branch cut of the complex plane before application of the Padé approximation. Denoting the rotation angle by , the representation of the square root is which, after expansion into a Padé series according to (6), yields where the complex Padé coefficients are given by Note that, in (11), is an approximation to one. This approximation gets better the more terms are used in the sum. However, for a finite number of terms , this approximation is always imperfect. Therefore, it is more practical to directly use . We will use this value for in the following derivations.
2.3. Fourier Finite Difference Migration
The complex Padé approximation (9) allows more stable implementations of not only FD but also FFD migration [20]. The derivation of 3D complex Padé FFD (CPFFD) migration is very similar to the original derivation of Ristow and Rühl [2]. It starts from the difference between the square root of (5) and a corresponding one where the velocity has been replaced by a constant reference velocity , namely, Expanding both square roots in (12) in complex Padé series according to (9), we find where we have used the notations
Joining the two series into one, expanding the fractions into Taylor series, and grouping the terms of equal power leads to Since this expression is, up to second order, equivalent to a Taylor series expansion of a Padé expression of the form where .
Using this approximation, the one-way wave equation (5) can be represented as which is the complex Padé equivalent to standard FFD migration [2].
As seen above, the theoretical value of obtained from this expansion is . However, 2D numerical experiments of Amazonas et al. [20] indicate that other expressions for can produce better results. For small contrasts, they suggest and, for high fidelity up to high propagation angles, they propose .
2.4. Implementation
To solve (18), we separate it into a set of differential equations. The first two terms provide the equations the analytic solutions of which are The remaining terms of (18) from the Padé series are represented by differential equations:
Discretizing these differential equations using a Crank-Nicolson FD scheme, we obtain where . Equation (22) means that the following implicit equation needs to be solved: We still need to discretize of the derivatives in the horizontal coordinates, that is, replace the differential operator by its difference operator where the matrices and represent difference operators for the second derivatives in and . For simplicity, we choose second-order difference operators, that is, for and with and denoting the number of grid points in the and directions. The resulting difference equation equivalent to differential equation (21) reads where is the identity matrix and is the matrix formed by the elements at a fixed depth level . Moreover, is the complex matrix with elements with , , and , and denotes the matrix obtained when replacing by in (27).
2.5. Two-Way Splitting
In the technique called two-way splitting, also known as alternating-directions-implicit (ADI) method [13, 24, 25], one substitutes (26) by its approximate factorized form This equation has the advantage of being solvable in two 2D steps. Under the assumption that the inverse operator in , , commutes with both operators in , equation (28) can be rewritten as where the intermediate value , defined as can be found solving the system The advantage of splitting is in the computational cost. While the numerical solution of (26) requires the solution of a system of size , the cascaded solution of (31) and (29) demands only the solution of systems of size , followed by systems of size . Since all these systems are tridiagonal, there are very efficient ways to solve them, which makes the splitting technique orders of magnitude faster than the solution of the original 3D system.
On the other hand, this procedure also has disadvantages. The biggest one is the introduction of numerical anisotropy into the propagation of the wavefield, because the numerical error increases with the azimuth between the propagation direction and the coordinate directions. This degrades the migrated image, introducing errors in the positioning of steeply dipping reflectors.
2.6. Splitting in More Directions
To overcome the problems with numerical anisotropy, Ristow and Rühl [5] proposed to generalize the technique to splitting into more than two directions. The idea is to approximate the 3D square-root operator by a sequence of 2D operators in different directions. In practice, most uses rely on three, four, or six directions to avoid symmetry problems. The unknown coefficients of these 2D operators are obtained from Taylor series expansions or by optimization techniques.
The multiway splitting form of the 2D Padé operators is given by the complex Padé expansion of the square-root operator in (5) for multiple directions: where is the number of directions. Moreover, are the derivative operators in the splitting directions, that is, , with (for , and ) being the azimuth of the rotated direction.
There are two ways of obtaining the unknown coefficients and in (32). One way is, as detailed above for the full 3D case, by Taylor series expansion of the fractions and comparison of the result with the direct Taylor series expansion of the square root. Alternatively, in practice, most often optimization techniques are employed to find optimum coefficients that minimize the numerical anisotropy for a certain range of medium velocities and a given reference velocity within the range of interest of propagation angles.
In conventional implementations of multiway splitting, operators (32) are applied in sequence at one single depth level before proceeding to the next one. In this paper, we apply the differential operators and for and , that is, the derivatives in the and directions, at one depth level, and leave the application of the operators and for and , that is, the derivatives in the diagonal directions, to the next depth level. In this way, we simulate four-way splitting, but with practically the same cost as conventional two-way splitting.
2.7. Stable FFD Migration and FFDPI Migration
We compare our results of CPFFD migration to another stable FFD migration technique, FFDPI migration [14]. It is based on a correction to the FFD method that avoids stability problems. To derive it, Biondi [14] starts from the real version of (18). In our notation, he rewrites the last part of the operator inside the summation as This representation of the operator corresponds to the differential equations This form of the differential equation can be implemented in a stable way by the realization of the product as a symmetrical matrix product , where is a diagonal matrix containing the values of the square root of . The remaining factor can also be represented as a product with a diagonal matrix.
In 3D, after two-way splitting, the resulting difference equation is approximated by the system where and are diagonal matrices with elements Moreover, matrices and are given by where and are diagonal matrices with elements and is given by According to Biondi [14], for the equivalence between the original differential equation and its stable discretization to be guaranteed, it is essential that is constant on the depth level under consideration, that is, that the reference velocity is larger or smaller than all model velocities at the current depth level.
Biondi [14] proves that this implementational correction stabilizes FFD migration, even in the presence of strong lateral velocity contrasts and for reference velocities larger than the medium velocity. In this way, this version of the FFD method possesses the necessary characteristics to be utilized as the main part of a precise and efficient high-angle wave-equation migration method. To attain a desired precision, one can interpolate between wavefields obtained for a sufficiently dense set of reference velocities. Theoretically, this allows to obtain arbitrary precision by increasing the number of reference velocities. The structure of this FFDPI method is similar to PSPI [19] and to the extended split-step method [26].
3. Numerical Experiments
3.1. Tests in a Homogeneous Medium
To study the numerical anisotropy of FFD migration operators after splitting, we calculated impulse responses for zero-offset migration in a homogeneous medium with velocity 2.5 km/s. The source pulse was a Ricker wavelet with central frequency 25 Hz, with its center positioned at an arrival time of 1.12 s. The migration grid was m, and m. All our examples used a complex Padé implementation of FFD migration with 3 terms in the series.
The top part of Figures 1(a)–1(d) show one vertical and three horizontal cuts through the reference impulse response, obtained with phase-shift migration using the true medium velocity. The amplitude decay at high propagation angles is caused by the source implementation, which did not use the amplitude correction of Wapenaar [27]. The red line in Figure 1(a) indicates the true theoretical position of the event, given by the half circle , where is the observation time of the event in the data, here 0.56 s, and is the source position, here the centre of the image, that is, m. The noncircular appearance of this line is due to the overstretched vertical axis. For a better comparison, we will present all other impulse responses below in the same way.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)
Figures 1(e)–1(h) show the corresponding four cuts of the impulse response of complex Padé FFD migration using conventional two-way splitting. Here, the value of the reference velocity was chosen as m/s, that is, . We observe the well-preserved circular shape of the impulse response in the deepest horizontal cut (Figure 1(i)), that is, for propagation directions close to the vertical axis. However, the shallow and, principally, medium horizontal cuts reveal a visible deformation, indicating the loss of quality for higher propagation angles. Also note the amplitude loss in the directions of the coordinate axes that are visible in the shallow and medium horizontal cuts. In the vertical cut (1(a)), only a slight deformation from circular shape is visible, which is due to the cut being within the coordinate plane, where the errors are the smallest. The amplitude loss for high propagation angles reflects the quality of the three-term Padé approximation. Note that the observed behaviour will be emphasized in media with strong lateral variations, where much smaller values of will occur.
Figures 1(i)–1(l) show the impulse response of FFD migration using conventional four-way splitting. The circular shape of the impulse response has been nicely recovered by the application of the two additional differential operators in the diagonal directions. Also, the amplitude loss in the coordinate directions is no longer visible. Note that this image has about twice the computational cost of the one in the center.
Figures 2(a)–2(d) show the impulse response of CPFFD migration using alternating four-way splitting, that is, two-way splitting in the coordinate directions at one depth level and in the diagonal directions at the next depth level. It is hard to spot any difference to the result of complete four-way splitting (Figures 1(i)–1(l)). The circular format of the operator is almost perfect, and even the slight amplitude loss along the coordinate axes is as well recovered as by complete four-way splitting. Note that this image has about the same computational cost as the one obtained with conventional two-way splitting (Figures 1(e)–1(h)).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)
Figures 2(e)–2(h) shows the impulse response of FFDPI migration using conventional two-way splitting, with interpolation between wavefields obtained for and . We chose these values to reflect the fact that, for FFDPI, generally reference velocities closer to the medium velocity are available for interpolation. We observe a good preservation of the circular shape, particularly in the horizontal cuts. In the vertical cut, we note that the wavefront lags slightly behind the true position, starting already at rather low propagation angles of about . The amplitude decay for high propagation angles is reduced as compared to FFD, probably because the reference velocities are closer to the medium velocity than in the previous examples. Finally, the shallowest cut exhibits some numerical dispersion, causing a distorted pulse shape.
Figures 2(i)–2(l) show the impulse response of FFDPI migration using alternating four-way splitting, with interpolation between wavefields obtained for and . Almost no improvement over the conventional two-way splitting result (center) is visible.
3.2. Tests in an Inhomogeneous Medium
For a more realistic test of the different splitting techniques for FFD migration, we calculated zero-offset impulse responses for the EAGE/SEG salt model. Here, we used a seismic pulse in the centre of the model, described by a Ricker wavelet with central frequency of 15 Hz, dislocated by s, and a migration grid with m. To avoid spurious events from the spike reflectors, we regularized the model using a median filter. This is necessary because the EAGE/SEG salt model uses artificially high velocity values to define the reflectors.
We represent the results by vertical cuts parallel to the - plane at km and km and parallel to the - plane at km and km, as well as horizontal cuts at depths m, km, km, and km. Figures 3 and 4 show these cuts through the EAGE/SEG salt model after filtering.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
Figures 5 and 6 show the corresponding cuts through the impulse response of FFD migration with two-way splitting and Figures 7 and 8 those of FFD migration with alternating four-way splitting. The differences between these sets of figures are due to numerical anisotropy, which is not always easy to see at this scale. The most visible difference is the one between Figures 6(a) and 8(a). The circular shape of three quarters of the wavefront is well preserved in Figure 8, while visibly distorted in Figure 6. Similar distortions are present in the other figure parts. Some events, particularly in the diagonal directions, are slightly more advanced in Figures 7 and 8 than in Figures 5 and 6. Also, some amplitude differences are visible. We refrain from presenting the results of complete four-way splitting, because they look virtually identical to those in Figures 7 and 8.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
For comparison, Figures 9 and 10 show the impulse response of FFDPI migration with two-way splitting. Since the theory of Biondi [14] is only formulated for a single term of the Padé series, so is our implementation. Because of the strong dependence of FFDPI on reference velocities not too far from the true model velocity, this numerical experiment needed 10 reference velocities. For being a very robust method, the impulse response is not subject to any instabilities, even with the reference velocities being still a bit far from the medium velocities. This remains true even for less reference velocities, though the image quality degrades considerably. Because of the need for a rather large number of reference velocities, FFDPI is a rather expensive method. In our implementation, it used about three times the computational time of alternating four-way FFD.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
Even for this experiment with 10 reference velocities, we still see some effects of numerical dispersion in Figures 9 and 10. Also, the results still exhibit quite visible differences to Figures 7 and 8. Since we have at this time no 3D reverse-time migration available, it is hard to tell which results are better positioned. Visual inspection and comparison to results of FD migration (not shown here) make us believe that the FFD results are more reliable than the FFDPI results with 10 reference velocities. More accurate results can be obtained by further increasing the number of reference velocities.
4. Conclusions
In this paper, we have implemented 3D versions of complex Padé FFD (CPFFD) and FFD plus interpolation (FFDPI), which have proven to be more stable in the presence of strong lateral velocity contrasts than other FFD migration implementations. For CPFFD migration, we have compared the effects of different ways of directional splitting and compared its results to those of FFDPI migration. Alternating four-way splitting, that is, applying the differential operators in the coordinate directions at one depth level and in the diagonal directions at the next depth level, proved to be an improvement over conventional two-way splitting at almost no extra cost. Although this procedure is theoretically less accurate than complete four-way splitting, that is, all four directions applied at all depth levels, our numerical results were of comparable quality. Extensions of the alternating splitting technique can be thought of like eight-way splitting where the remaining directions are covered two by two in the next two depth steps.
From our numerical tests with splitting the CPFFD and FFDPI migration operators, we conclude that FFDPI migration is the most robust of the tested methods. Even implemented only using two-way splitting, it did show only a fair amount of numerical dispersion, but no visible numerical anisotropy. However, for practical use, FFDPI is a rather expensive method because it needs a large number of reference velocities to function with acceptable precision. Thus, for a more economic migration with acceptable image quality, alternating four-way splitting in FFD migration is an interesting alternative.
One minor problem of multiway splitting should be mentioned. The differential operator in the diagonal directions can cause aliasing effects because of the fact that the grid spacing in this direction is by a factor of larger than in the coordinate directions. Off-diagonal directions may complicate things further, because they require resampling.
Acknowledgments
This work was kindly supported by the Brazilian agencies CAPES, FINEP, and CNPq, as well as Petrobras and the sponsors of the Wave Inversion Technology (WIT) Consortium.