# Three-Dimensional Complex Padé FD Migration: Splitting and Corrections

Authors: D. Mondini; J. C. Costa; J. Schleicher; A. Novais
Journal: International Journal of Geophysics
 (2012)

---

## 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.


---

## Body

## 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 x and y 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 timet=0. 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](1)∇2p(x,t)=∂2p(x,t)∂x2+∂2p(x,t)∂y2+∂2p(x,t)∂z2=1v2∂2p(x,t)∂t2,
where ∇2 is the Laplacian operator, p(x,t) is the scalar wave field, x=(x,y,z), and v is the speed of the medium, considered constant. We define the Fourier transform in horizontal coordinates x and y and time t as(2)P(kx,ky,z,ω)=∫-∞∞p(x,t)ei(kxx+kyy-ωt)dtdxdy,
and its inverse in the same variables as(3)p(x,t)=1(2π)3∫-∞∞P(kx,ky,z,ω)e-i(kxx+kyy-ωt)dωdkxdky.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 x is well described by a constant speed v, but the value of this constant velocity depends on x, 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(4)[∂∂z-iωv(x)1-v(x)2ω2(kx2+ky2)]P(kx,ky,z,ω)=0.In the space domain, this equation can be represented formally as(5)∂P(x,ω)∂z=iωv(x)1+v(x)2ω2(∂2∂x2+∂2∂y2)P(x,ω),
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 integersL and M denoting the degrees of the numerator and denominator polinomials, and are represented by the notation [L/M].  Explicitly, an Padé approximant is defined by(6)[LM]=PL(x)QM(x),L,M≥0.Using the method of Padé approximants in (5), we can rewrite the square root operator as(7)1+Z≈1+∑n=1NanZ1+bnZ,
where the coefficients an and bn are real Padé coefficients [5] given by(8)an=22N+1(sin⁡2nπ2N+1),bn=cos⁡2(nπ2N+1).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 90°, the argument Z becomes very close to −1. And when Z<-1, 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(9)1+Z≈C0+∑n=1NAnZ1+BnZ,
where the complex Padé coefficients are given by(10)An=ane-iα/2[1+bn(e-iα-1)]2,Bn=bne-iα1+bn(e-iα-1),C0=eiα/2[1+∑n=1Nan(e-iα-1)[1+bn(e-iα-1)]]≈1.
The constant C0 is an approximation to 1 that is the better the more terms N are used. Since it does not depend on the argument Z, 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 an and bn by the corresponding complex coefficients An and Bn, resulting in(11)∂P(x,ω)∂z=iωv(x)[1+∑n=1NAn(v2(x)/ω2)(∂2/∂x2+∂2/∂y2)1+Bn(v2(x)/ω2)(∂2/∂x2+∂2/∂y2)]×P(x,ω).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(12)∂P(x,ω)∂z=iωv(x)P(x,ω),
followed by equations resulting from the terms inside the summation, that is,(13)∂P(x,ω)∂z=iωv(x)[An(v2(x)/ω2)(∂2/∂x2+∂2/∂y2)1+Bn(v2(x)/ω2)(∂2/∂x2+∂2/∂y2)]P(x,ω),
for n=1,…,N.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 v(x)=v(x,y) within a layer of size Δz. Then, we can express the wavefield at level z+Δz depending on the wavefield at the previous level, z, as(14)P(x,y,z+Δz,ω)=P(x,y,z,ω)exp⁡[iωv(x)Δz].The differential equations associated with the summation are given by(15)∂P(x,ω)∂z+Bnv2(x)ω2(∂2∂x2+∂2∂y2)∂P(x,ω)∂z=-Anv(x)iω(∂2∂x2+∂2∂y2)P(x,ω).
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 x or y 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 velocityv(x) by a reference velocity vr 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Δx=Δy=Δz=20m (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 t=1.1 s in time and x=6.75km, y=6.75km, and z=0km in space.To represent the results, we use vertical cuts parallel to theyz and xz planes at x=y=8.32km and horizontal cuts at z=1.1km and z=2.7km 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 (x,y) 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.Two horizontal and two vertical cross-sections across the EAGE/SEG salt model atx=8.32km (a), y=8.32km (b), z=1.1km (c) and z=2.7km (d).
(a)(b)(c)(d)Unless otherwise specified, we use in our numerical tests the following parameters: three terms (N=3) in the complex Padé approximation, rotation angle α=45°, 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.Impulse response obtained by FFDPI migration represented by 4 cross-sections atx=8.32km (a), y=8.32km (b), z=1.1km (c), and z=2.7km (d).
(a)(b)(c)(d)In Figure3 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).Impulse response obtained by FD migration using the complex Padé approximation without Li correction (cross-sections as before).
(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.Figures4 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.Impulse response obtained by FD migration using the real Padé approximation with 1 term and Li correction (cross-sections as before).
(a)(b)(c)(d)Impulse response obtained by FD migration using the complex Padé approximation with 1 term with rotation angle of15° and Li correction (cross-sections as before).
(a)(b)(c)(d)Figures6 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.Impulse response obtained by FD migration using the real Padé approximation with 2 terms and Li correction (cross-sections as before).
(a)(b)(c)(d)Impulse response obtained by FD migration using the complex Padé approximation with 2 terms with rotation angle of25° and Li correction (cross-sections as before).
(a)(b)(c)(d)In Figures8 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.Impulse response obtained by FD migration using the real Padé approximation with 3 terms and Li correction (cross-sections as before).
(a)(b)(c)(d)Impulse response obtained by FD migration using the complex Padé approximation with 3 terms with rotation angle of45° and Li correction (cross-sections as before).
(a)(b)(c)(d)Figure10 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.Impulse response obtained by FD migration with rotation angle of45°, 3 terms, mean slowness, and Li correction (cross-sections as before).
(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. Figure11 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.Impulse response obtained by FD migration using alternating four-way splitting and without Li correction (cross-sections as before).
(a)(b)(c)(d)Impulse response obtained by FD migration using alternating four-way splitting and with Li correction (cross-sections as before).
(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 factor2 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.



---

*Source: 479492-2012-05-16.xml*