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. To reduce errors, the Li correction is applied at regular multiples of depth extrapolation increment. We compare the performance of splitting techniques in wave propagation for 3D finite-difference (FD) migration in terms of image quality and computational cost. We study the behaviour of the complex Padé approximation in combination with two- and alternating four-way splitting, that is, splitting into the coordinate directions at one depth and the diagonal directions at the next depth level. We also extend the Li correction for use with the complex Padé expansion and diagonal directions. From numerical examples in inhomogeneous media, we conclude that alternate four-way splitting is the most cost-effective strategy to reduce numerical anisotropy in complex Padé 3D FD migration.
1. Introduction
In three dimensions, migration methods based on solving the one-way wave equation, besides facing problems to image dip reflectors and handle evanescent waves, are still computationally expensive. For the problems of imaging dip reflectors and evanescent waves, we use the complex Padé approximation. Because the resolution of three-dimensional problem is computationally expensive, over the years various techniques have been developed in order to reduce costs and still maintain the quality of the migration method in use. A commonly used technique is splitting.
For the case of splitting in two directions, we face the problem of numerical anisotropy, that is, the migration operator acts differently in different directions, resulting in positioning errors of reflectors in the situation where the direction of the dip reflector is far from the directions of the migration planes. To correct this problem it is common to use the correction of Li [1]. Without changing the basic principle of applying subsequent 2D FD migrations in the and directions, the Li correction is an extrapolation of the residual field by a phase shift. In our extension using complex Padé, the modified Li correction also includes compensation for evanescent energy propagation. When splitting is applied alternately in four directions (the horizontal coordinates and the diagonals), we may still face problems of numerical anisotropy and, consequently, of positioning errors of steeply dipping reflectors. Therefore, we also tested the application of a Li correction adapted to this case.
Our goal in this work is to evaluate the behaviour of 3D FD migration operators using the complex Padé approximation. We investigate the technique of splitting into two or four alternating directions as well as the corresponding Li corrections in order to determine an effective strategy to reduce computational cost and numerical anisotropy. For that purpose, we compare the impulse responses obtained by FD migration in the SEG/EAGE salt model.
2. Theory
According to the migration imaging condition, a migration consists only in repositioning the seismic wave recorded at the source position to depth at time . Thus, we see that we need only to propagate the waves downward from the sources and the receivers into the earth to perform a migration. Therefore, it is useful to separate the wave equation into two one-way wave equations that describe upward and downward propagation separately [2].
The acoustic wave equation for a homogeneous medium is given by [3] where is the Laplacian operator, is the scalar wave field, , and is the speed of the medium, considered constant. We define the Fourier transform in horizontal coordinates and and time as and its inverse in the same variables as
Applying the Fourier transform as defined in (2) to (1) and considering that the velocity is locally constant, that is, the propagation of waves at each location is well described by a constant speed , but the value of this constant velocity depends on , we can separate the resulting equation into two equations describing the up and downgoing wavefields, respectively. The one-way wave equation for downgoing waves reads
In the space domain, this equation can be represented formally as which is the downgoing one-way wave equation for media with moderate-velocity variations.
3. Complex Padé Approximation
To actually use (5) in migration, we need to approximate the square-root operator by some numerically executable expression. Using a Taylor approximation may present some difficulties inherent in such series. Often the convergence of the series is extremely slow, or else, its radius of convergence does not include areas of particular interest to the problem being studied. On the other hand, Padé series enable us to start from a power series and obtain much more information than the series itself can provide us directly [4].
Padé approximants are rational functions, that is, quotients of two polynomials, which represent an expansion. These approximants are characterized by two positive integers and denoting the degrees of the numerator and denominator polinomials, and are represented by the notation . Explicitly, an Padé approximant is defined by
Using the method of Padé approximants in (5), we can rewrite the square root operator as where the coefficients and are real Padé coefficients [5] given by
The Padé approximation [6] is widely used in practice, because it allows an efficient numerical implementation and a larger error order than other approximation methods. However, when looking for images with good accuracy for wide angles of propagation, (7) becomes inappropriate. The reason is that for angles near , the argument becomes very close to −1. And when , the left side of (7) is complex while the right side is still real, causing an inconsistency in the approximation [7]. Physically, this means that the representation of (7) cannot properly handle evanescent waves. This feature is responsible for the unstable behaviour of the algorithm in the presence of strong lateral velocity variations.
In the complex plane, a Padé expansion with real coefficients corresponds to approximations of the square root with branch cut along the negative real axis. To overcome the problems with evanescent waves, Milinazzo et al. [8] proposed the introduction of a rotation angle , which rotates the branch cut, improving the representation of the evanescent part of spectrum and therefore the stability of the approach.
Considering the rotation of the branch cut on the complex plane, the representation of the square root has the form where the complex Padé coefficients are given by The constant is an approximation to 1 that is the better the more terms are used. Since it does not depend on the argument , it can be immediately replaced by 1 in the algorithm.
4. FD Migration with the Complex Padé Approximation
This method approximates the operator of the one-way wave equation (5) by a Padé series [5, 6]. Its complex version is obtained replacing the real coefficients and by the corresponding complex coefficients and , resulting in
Like its real version, this operator can be implemented iteratively, considering each term of the summation independently. Thus, we first need to solve the differential equation followed by equations resulting from the terms inside the summation, that is, for .
To solve this system of equations, we consider each solution of the previous equation the initial condition for the following equation. The first equation (12) has an analytic solution for models for which the velocity satisfies within a layer of size . Then, we can express the wavefield at level depending on the wavefield at the previous level, , as
The differential equations associated with the summation are given by The advantages of using the complex version are that this technique does not have limitations towards the velocity model variations. It can handle evanescent waves more appropriately. The downside remains the same as in the real version of the FD method, which is the difficulty of imaging steep-dip reflectors.
5. Splitting and Li Correction
The numerical solution of the differential equation (15) involves the inversion of a huge matrix, which is a very costly process. Individual solution in the coordinate directions [suppressing the or derivatives in (15)] is orders of magnitude less expensive. So, the idea arose to approximately execute the full 3D operator of (15) as a sequence of 2D processes [9]. This approach is called directional splitting.
Generalizing this idea, Ristow and Rühl [10] proposed splitting in various directions, which consists in approximating a 3D operator in a sequence of 2D operators, usually along the lines of horizontal coordinates, the inline and crossline directions and some additional directions away from the coordinate axis. The price to pay for the increased efficiency of the method is reduced accuracy in the directions not used for the 2D operators.
6. Correction for Two-Way Splitting (Li Correction)
To compensate the errors caused by the use of splitting into the coordinate directions (below referred to as two-way splitting) and still preserve the efficiency of the FD method, Li [1] proposed the application of a phase correction operator, implemented either as a phase-shift term or using the PSPI method. This difference operator is obtained by evaluating the error between the true and approximate operators.
The idea of this correction is to carry out the conventional 2D FD migration in both directions where the splitting was done plus a residual wavefield extrapolation by the phase-shift method (when the lateral velocity variation is small) or the PSPI method (when the lateral variations of the wave speed are more significant). The modified Li correction for complex Padé FD also includes compensation for evanescent energy propagation.
Because the error is supposedly small in each step of the extrapolation, the effect of the compensation process is similar to the residual migration, in which case the waves propagate very little at each step. Therefore, it is reasonable to replace the true velocity by a reference velocity and apply the correction only at a reduced number of depth steps. To compute representative reference velocities we need to average the velocity field along the depth interval between successive applications of the Li correction. We investigated the most effective way how to compute this representative velocity field by averaging the slowness, the velocity, or the cube of the velocity, as proposed by Li [1]. We also investigated the lateral smoothing of this average field to reduce strong lateral variations in the Li correction.
The average velocity field is the input to compute reference velocities for Li correction. For the choice of these reference velocities, several methods exist. In this work, we employ a modified version of the method of Lloyd [11]. The main idea behind this method is to iteratively improve the values of a chosen set of reference velocities by looking at statistical values in each region (mean, median, and variance) and then change the region boundaries in each iteration as an attempt to find the best solution based on some criterion. In our modified version, postprocessing is applied to the Lloyd velocities at the end of the process. Of every three reference velocities within an interval of 0.3 km/s, the central one is eliminated. Moreover, at large velocity gaps in the model, the Lloyd reference velocities tend to fall into the center of the gap. We replace these values by the two model velocities that border the gap.
7. Results
We tested these implementations of wave propagation for FD migration by their impulse responses in the EAGE/SEG salt model. We used the following modelling parameters: The source wavelet is a Ricker pulse with central frequency of 15 Hz; the migration grid spacing is m (to avoid numerical dispersion, we treat the model as if the spacing was 10 m); the delta pulse to be migrated is symmetrically positioned with its centre at s in time and km, km, and km in space.
To represent the results, we use vertical cuts parallel to the and planes at km and horizontal cuts at km and km depths (see Figure 1). In order to use Lloyd’s method for choosing the reference velocities for each level, we calculate the average slowness using all points of the velocity model. Because it is an iterative method, the number of reference velocities is variable, limited to at most 10 velocities at each level.
(a)
(b)
(c)
(d)
Unless otherwise specified, we use in our numerical tests the following parameters: three terms () in the complex Padé approximation, rotation angle , and application of the Li correction at every 6 steps. We compare the results varying these parameters individually. As a reference, we use the result of the FFDPI method [12] to evaluate our results (see Figure 2). Despite the numerical dispersion visible in these figures, they show no signs of numerical anisotropy, so that the positioning of the events can be considered reliable.
(a)
(b)
(c)
(d)
In Figure 3 we see the result of complex Padé FD migration with three terms and without Li correction. In the horizontal sections, we see a deviation from circular to diamond shape and the appearance of artifacts, indicating the strong numerical anisotropy caused by splitting. We also note that the wavefronts lag behind their true position (compare to Figure 2).
(a)
(b)
(c)
(d)
Next, we investigate the dependence on the number of terms in real and complex FD migration. As a general observation, we find that the optimum rotation angle depends on the number of terms used. Below, we show only the results with the best rotation angle.
In another test, we investigated different ways to average the velocity field in the depth interval between two applications of the Li correction. We tested average velocity, slowness, and the cube of velocity. The best results were produced when averaging the slowness. We also evaluated the effect of laterally smoothing this average velocity field with a moving average filter. This smoothing it did not improve the Li correction. The applications of Li correction below use average slowness to estimate the representative velocity field in the depth interval between two corrections.
Figures 4 and 5 show the impulse responses of real and complex Padé FD migration, both with 1 term and Li correction. Comparing Figures 4 and 5 we note the appearance of artifacts in the real version in the vertical cross-sections and blurred events in the complex version. The artifacts present in the shallowest horizontal cut are a consequence of applying the Li correction. However, there are less artifacts in the complex version than in the real one. The comparison with the results in Figure 3 shows that Li correction actually improves the positioning of the events. Delays in the diagonal directions were fixed, recovering an almost perfect circular shape.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
Figures 6 and 7 show the corresponding results with 2 terms and Li correction. Comparing Figures 6 and 7 with Figures 4 and 5 we note that with 2 terms we have less artifacts in the real version and more marked events in the complex version in the vertical cuts.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
In Figures 8 and 9 we see the corresponding results obtained with 3 terms. Comparing Figures 4 and 5 we note that we have less artifacts than when we use 1 term in the approximation. The 3-term real result attains approximately the same quality as the 2-term complex result.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
Figure 10 depicts the results when using the mean velocity instead of the mean slowness for the reference velocity computation. As expected, the result is of inferior quality to that obtained with the mean slowness (Figure 9), because the latter corresponds to the real representation we can see in the wave equation.
(a)
(b)
(c)
(d)
In addition to the above experiments with two-way splitting, we also tested alternating four-way splitting, where the splitting is carried out in the coordinate directions at one level and in the diagonal directions at the next level. Apart from minor effects, alternating four-way splitting has the same computational cost as two-way splitting. Figure 11 shows the results. We see that the problem of numerical anisotropy has improved but not completely solved. Although we do not see artifacts in the horizontal cuts, we see a deviation from circular to octagonal shape. Other remaining positioning errors are evident when comparing Figure 11 with the corresponding Figure 2 obtained by FFDPI method. Differently from what has been reported for FFD migration by Costa et al. [13], alternating four-way splitting without Li correction is clearly not sufficient for FD migration even for moderately varying velocity. To fix the remaining anisotropy we also applied the Li correction (see Figure 12). The result has improved considerably.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
8. Conclusions
From the numerical tests, we conclude that FD migration with two-way splitting plus Li correction showed an excellent recovery of circular shape but still showed some remaining phase errors. As a drawback, it must be noted that as a consequence of the Li correction, some migration artifacts emerged or were reinforced.
When comparing the real and complex versions of the Padé approximation, we found few differences. The use of fewer terms in the Padé expansion affected the real version more strongly than its complex counterpart. With respect to the rotation angle of the branch cut in the complex Padé approximation, we found that the fewer number of terms are used the smaller the rotation angle must be.
The application of splitting in more than two directions is useful to reduce the numerical anisotropy but is not cost effective. Moreover, in the diagonal directions aliasing effects may arise because the distances of the grid points are a factor higher. Other directions outside of the axes and diagonals complicate the problem because they demand resampling. One way to reduce the cost of splitting into four directions is the application in alternating directions, the coordinate directions at one depth step, and the diagonals at the next. This alternating application has the same cost of splitting in two directions but has a very high quality. In our experiments, the quality was equivalent to full four-way splitting. Thus, the numerical results in four directions presented in this paper exclusively use this alternate implementation. For not too deep reflectors, the remaining numerical anisotropy is acceptable. For larger migration distances, regardless of the splitting into two-way or alternating four-way directions, Li correction should not be omitted if accurate positioning is desired.
We also investigated different strategies for the application of the Li correction. Best results were achieved when using the average slowness in the depth interval between successive applications, field interpolation using velocity, and reference velocities chosen by Lloyd’s method.
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.