TY - JOUR
AU1 - Costa, Carlos A, N
AU2 - Campos, Itamara, S
AU3 - Costa, Jessé, C
AU4 - Neto, Francisco, A
AU5 - Schleicher,, Jörg
AU6 - Novais,, Amélia
AB - Abstract Conventional implementations of 3D finite-difference (FD) migration use splitting techniques to accelerate performance and save computational cost. However, such techniques are plagued with numerical anisotropy that jeopardises the correct positioning of dipping reflectors in the directions not used for the operator splitting. We implement 3D downward continuation FD migration without splitting using a complex Padé approximation. In this way, the numerical anisotropy is eliminated at the expense of a computationally more intensive solution of a large-band linear system. We compare the performance of the iterative stabilized biconjugate gradient (BICGSTAB) and that of the multifrontal massively parallel direct solver (MUMPS). It turns out that the use of the complex Padé approximation not only stabilizes the solution, but also acts as an effective preconditioner for the BICGSTAB algorithm, reducing the number of iterations as compared to the implementation using the real Padé expansion. As a consequence, the iterative BICGSTAB method is more efficient than the direct MUMPS method when solving a single term in the Padé expansion. The results of both algorithms, here evaluated by computing the migration impulse response in the SEG/EAGE salt model, are of comparable quality. 93.85.Rt, 93.85.Bc 3D seismic imaging, wave-equation, sparse-solvers 1. Introduction Efficient implementations of 3D finite-difference (FD) use the splitting technique along horizontal coordinates in order to avoid solving large-band linear systems (Claerbout 1985) at each downward continuation step. Splitting the migration operator along inline and crossline directions reduces 3D migration to a sequence of 2D downward continuation steps, thereby making the algorithm highly efficient. However, splitting introduces numerical anisotropy, which can damage the image of reflectors (Brown 1983). Several strategies are used to remedy this problem, usually by assuming a homogeneous medium and using wavefield interpolation to handle lateral velocity variations (Li 1991). Other ideas to mitigate the effects of splitting are to increase the number of splitting directions (Ristow and Rühl 1997) or to alternate splitting directions and add an interpolation step (Wang 2001). However, such correction methods add again to the cost and introduce additional errors and artefacts. Exploring another idea to avoid the expensive direct solution of these large systems, Cole (1989) investigated the use of iterative methods to implement FD migration without splitting and evaluated the performance of the overrelaxation, Jacobi, and Gauss–Seidel methods (Iserles 1996). He reported the poor conditioning of the linear system matrix for low frequencies, which degraded the performance of iterative solvers and led to a prohibitively high computational cost. Along a similar line, Nichols (1997) investigated the performance of the conjugate gradient (CG) method and its dependence on the choice of the initial value for preconditioning. He found a degradation of performance of the CG method for low frequencies but reported that the performance improved when the propagation of energy associated with high angles and evanescent modes was attenuated. In this respect, it is interesting to note that Amazonas et al (2007) improved the quality of migrated images by representing the downward continuation operator using the complex Padé approximation (Millinazzo et al1997), which is essentially a means of attenuating the propagation of evanescent wave modes. Moreover, new iterative methods to solve complex linear systems are available that allow the solution of very large linear systems to be computed. For example, the stabilized biconjugate gradient (BICGSTAB) method can solve poorly conditioned complex linear systems more efficiently than the CG method (van der Vorst 1992). Other approaches to improve conditioning include helical preconditioning (Rickett et al1998, Zhang et al2000), and domain decomposition (Fei and Etgen 2002). An alternative to the iterative approach are direct solvers. Here, we specifically use the multifrontal massively parallel direct solver (MUMPS), which can efficiently solve complex large-band linear systems using Gaussian decomposition (Amestoy et al2001, 2006). Costa et al (2011) and Mondini et al (2012) demonstrated the benefits of using the complex Padé approximation for operator splitting in FD and Fourier FD migration. This approximation greatly reduces migration artefacts from instabilities due to strong velocity contrasts. Costa et al (2011) showed that in moderately complex media, FFD migration with four-way splitting using the complex Padé approximation can produce an acceptable image quality without the need for any further correction. On the other hand, FD migration, even with four-way splitting and the complex Padé approximation, always requires Li correction (Mondini et al2012). In this paper, we investigate whether the limitations observed by Nichols (1997) can be overcome when using the complex Padé approximation in the implementation of 3D downward continuation by FD without splitting. Our numerical experiments in homogeneous media and in the SEG/EAGE salt model (Aminzadeh et al1995) indicate that the complex Padé approximation is an effective preconditioner for 3D FD migration, so that the iterative BICGSTAB approach outperforms the massively parallel MUMPS solver when using a single term of the Padé expansion. If more terms are needed to better model high dip events, the iterative method fails. 2. Theory 2.1. Implicit finite-difference downward continuation The one-way wave equation for downward continuation in the space-frequency domain is ∂P(x,ω)∂z=-iωc(x)1+DP(x,ω),1 where ω is the angular frequency and x = (x, y, z) is the position with x and y denoting the coordinates along inline and crossline directions and z denoting depth. Moreover, P(x, ω) is the wavefield spectrum and c(x) represents the wave speed. Finally, D represents the differential operator D=c2(x)ω2∂2∂x2+∂2∂y2.2 The square-root operator in (1) can be expanded by using the complex Padé series (Millinazzo et al1997, Amazonas et al2007): ∂P∂z=-iωc(x)C0+∑n=1NAnD1+BnDP(x,ω),3 where N indicates the order of the Padé expansion and C0, An, and Bn are the complex coefficients C0=eiθ/21+∑n=1Nan(e-iθ-1)1+bn(e-iθ-1),An=ane-iθ/2[1+bn(e-iθ-1)]2,Bn=bne-iθ1+bn(e-iθ-1).4 Here, θ indicates the angle of rotation for the branch cut of the square root in the complex plane, and an and bn are the real Padé expansion coefficients (Bamberger et al1988) an=22N+1sin2nπ2N+1andbn=cos2nπ2N+1.5 The constant C0 does not depend on D and is actually an approximation to one that gets better with an increasing number of terms N used in the expansion. Therefore, we can set C0 = 1. Note that the complex Padé expansion of Amazonas et al (2007) used here differs from the one proposed by Zhang et al (2003) (see also Zhang et al2004, 2005, 2007). In those works, the expansion is applied to the complete exponential term in the solution to the one-way wave equation (1), while in Amazonas et al (2007) and in the present work, the square root in (1) itself is expanded. Therefore, the restrictions that apply to the expansion given by Zhang at his coworkers do not apply here. As shown in Amazonas et al (2007), the square-root expansion is more stable and helps to reduce migration artefacts. The numerical solution of (3) is obtained by solving a sequence of differential equations (Claerbout 1985). At each continuation step Δz, the solution is computed sequentially for each term of the expansion, using the solution of the previous step as the initial condition for the next one. The first equation to be solved is ∂P0(x,ω)∂z=-iωc(x)P0(x,ω),6 followed by the N equations ∂Pn(x,ω)∂z=-iωc(x)AnD1+BnDPn(x,ω)(n=1,...,N),7 where the solution Pk(x, ω) (k = 0, …, N - 1) to each equation becomes the initial condition for the subsequent equation. This algorithm converges for small depth steps Δz. Using the Crank–Nicolson FD scheme (Iserles 1996) to approximate the derivative operators, equations (7) are approximated by the linear systems (I+α+D)Pk+1=(I+α-D)Pk,8 where α±≡Bn±iωΔzcAn,9 and D≡c2ω2Dx(Δx)2+Dy(Δy)2.10 Here, I is the identity matrix, Dx and Dy are the FD matrices representing the second derivatives along the x and y coordinates, Δx and Δy are the step sizes in these directions, and Pk and Pk + 1 represent the wavefield at the grid at successive depth levels. The matrix D has the dimension Nx × Ny. Though sparse, this matrix can have a large band. This increases the computation cost of solving system (8) using direct methods. An efficient and frequently used method to reduce this cost is splitting along inline and crossline directions, I+α+cωΔx2DxI+α+cωΔy2DyPk+1=I+α-cωΔx2DxI+α-cωΔy2DyPk.11 The solution of this linear system is a sequence of tridiagonal linear systems of dimensions Nx and Ny. Unfortunately, this approximation introduces numerical anisotropy through the incorrect mixed-derivative term. To reduce numerical anisotropy, Li (1991) proposed a correction in the frequency-wavenumber domain which requires wavefield interpolation in laterally heterogeneous media and Ristow and Rühl (1997) generalized the above splitting approach to allow for additional directions, the so-called multiway splitting. The corrections for splitting improve the image quality and effectively reduce numerical anisotropy. However, by their approximative nature, they do not fully remove all splitting effects and, in addition, even introduce further propagation artefacts. Moreover, they add to the computational cost. For these reasons, we are interested in alternative methods to solve system (8) that do not make use of splitting techniques. 2.2. Iterative methods The sparseness of linear system (8) makes it suitable for iterative methods (Iserles 1996). Cole (1989) investigated the performance of some iterative methods to solve (8) using the first-order real Padé expansion. He noticed the rather slow convergence, or even divergence, of methods like the overrelaxation, Jacobi, and Gauss–Seidel ones, when applied for frequencies below ω2 = 8c2bn/(Δx)2. Nichols (1997) used the (CG) method to solve (8) because of its better theoretical convergence properties, but also reported the poor performance of this method for low frequencies. However, he observed some improvement in the convergence speed of CG when evanescent modes were filtered out. Unfortunately, the performance degraded rapidly for ill conditioned matrices, because the CG method solves the normal equations. Since the days of Cole (1989) and Nichols (1997), iterative inversion techniques have improved. We investigate whether the limitations of conventional iterative methods as reported in the cited works still hold for more modern techniques. Taking into account the observation of Nichols (1997) that the performance of the CG method improves when the propagation of energy associated with high angles and evanescent modes is attenuated, we combine the complex Padé expansion (Amazonas et al2007) with the BICGSTAB algorithm (van der Vorst 1992) to implement downward continuation without splitting. We chose BICGSTAB because of its better convergence properties than the conventional CG method for complex matrices. 2.2.1. Conditioning The performance of iterative methods improves when the matrix to be inverted is strictly diagonally dominant (Iserles 1996). For a matrix A this means |Ajj|≥∑k=1NA|Ajk|j≠k(j=1,...,NA),12 where NA denotes the dimension of matrix A. For a homogeneous medium, we can verify when the matrix on the left side of (8) satisfies this condition. Its diagonal entries are Ajj=ωΔxc2-4Bn+iωΔx2cAn,13 and the only four nonzero off-diagonal entries in each row are all equal to Ajk=Bn+iωΔx2cAn(j≠k).14 Thus, this matrix is diagonally dominant if ω≥ωL=2cΔx[- Im {An}+( Im {An})2+2 Re {Bn}],15 where Im{An} and Re{Bn} represent the imaginary and real parts of An and Bn, respectively. This inequality determines the minimum frequency for which the matrix is strictly diagonally dominant. For frequencies above this limit, iterative methods should perform better than for lower frequencies. Although this result is only a sufficient condition, it sheds some light on the poor performance of iterative methods for downward continuation, particularly regarding their low-frequency behaviour. The limit (15) depends on the ratio between propagation velocity and grid interval. When this ratio increases, i.e., for small grid sizes, the convergence of iterative methods is slower. Another consequence from (15) regards the difference between the real and complex Padé expansions. In the real version, Im{An} = Im{Bn} = 0, and (15) reduces to ω≥ωLr=2cΔx2Bn.16 This implies that using the complex Padé expansion will improve convergence, and it will improve further the larger Im{An} is relative to Re{Bn}. This fact is corroborated in the top left part of figure 1, which shows how the normalized limit frequency ωL/ωLr decreases as a function of the branch-cut rotation angle θ (green line), for N = 1. The values at θ = 0° (real Padé) and θ = 45° are singled out by a blue and a red line, respectively. Figure 1. Open in new tabDownload slide Top left: normalized limit frequency, fL, versus rotation angle in degrees. Top right: number of iterations of BICGSTAB for each frequency for c/Δx = 75 Hz. Bottom left: iterations versus frequency for c/Δx = 150 Hz, Bottom right: iterations versus frequency for c/Δx = 450 Hz. Figure 1. Open in new tabDownload slide Top left: normalized limit frequency, fL, versus rotation angle in degrees. Top right: number of iterations of BICGSTAB for each frequency for c/Δx = 75 Hz. Bottom left: iterations versus frequency for c/Δx = 150 Hz, Bottom right: iterations versus frequency for c/Δx = 450 Hz. To evaluate how the actual behaviour of BICGSTAB relates to inequality (15), the remaining graphs in figure 1 display the number of iterations as a function of frequency for increasing values of the ratio c/Δx for real and for complex Padé coefficients with θ = 25° and θ = 45°. These graphs show that the complex Padé approximation markedly improves the convergence of BICGSTAB, particularly for large values of the ratio c/Δx. Increasing the branch-cut rotation angle from 0° (real approximation) to 25° strongly impacts the convergence, while the increase from 25° to 45°, though further reducing the number of iterations, has a less pronounced effect. 3. Numerical results 3.1. Homogeneous medium After the evaluation of the performance, we are interested in the quality of the resulting wavefield. As a first test, we computed the impulse response of the 3D FD migration operator in a homogeneous medium with and without inline and crossline splitting. The propagation velocity is 1500 m s-1. We used a single term of the Padé expansion (the complex versions with rotation angle θ = 25°) to approximate the migration operator. The impulse responses without splitting were computed using MUMPS and BICGSTAB. The grid dimensions are 676 × 676 × 210 along coordinates x, y, and z, respectively. The grid spacing is uniform and equal to 10 m. The volume injection source is located one grid spacing below the surface in the centre of the grid. The source signature is a Ricker pulse with peak frequency 25 Hz, and the sample rate is 8 ms. The resulting cross-sections of the 3D impulse responses at depth z = 1050 m and diagonal cross-sections are depicted in figures 2 and 3. The red circles in these images indicate the correct position. Figures 2 and 3 compare horizontal and diagonal sections through the wavefield for a homogeneous medium, computed with (a) splitting along the inline and crossline directions according to (11), (b) four-way splitting along the inline, crossline, and diagonal directions according to Ristow and Rühl (1997), and (c) four-way splitting with the correction of Li (1991) to the ones obtained without splitting using (d) real MUMPS, (e) complex BICGSTAB and (f) complex MUMPS, (8), at the same depth. The red circles in these images indicate the correct position. The noncircular format of the result of splitting in figure 2 is the consequence of numerical anisotropy because of the selected preferential directions. Moreover, we see that operator splitting also introduces noncausal artefacts that reach far beyond the correct wavefront. This is even true in the diagonal section obtained using four-way splitting (figure 3(c)), which is in the operator direction and should thus be exact. The high-frequency scatter in parts (a), (b), and (d) of figure 3 appears because of instabilities in the FD calculations. These artefacts are greatly reduced by the complex Padé approximation (parts (e) and (f) of figures 2 and 3). The improvement of the impulse response computed by MUMPS or BICGSTAB over the one computed using operator splitting is evident. The correct positioning and the absence of numerical anisotropy indicate the good performance of the BICGSTAB algorithm. Figure 2. Open in new tabDownload slide Horizontal section through the FD impulse response in a homogeneous medium. (a) Splitting along the inline and crossline directions; (b) four-way splitting in inline, crossline and diagonal directions; (c) four-way splitting plus Li correction; (d) real MUMPS; (e) complex BICGSTAB; (f) complex MUMPS. Figure 2. Open in new tabDownload slide Horizontal section through the FD impulse response in a homogeneous medium. (a) Splitting along the inline and crossline directions; (b) four-way splitting in inline, crossline and diagonal directions; (c) four-way splitting plus Li correction; (d) real MUMPS; (e) complex BICGSTAB; (f) complex MUMPS. Figure 3. Open in new tabDownload slide Diagonal section through the FD impulse response in a homogeneous medium. (a) Splitting along the inline and crossline directions; (b) four-way splitting in the inline, crossline, and diagonal directions; (c) four-way splitting plus Li correction; (d) real MUMPS; (e) complex BICGSTAB; (f) complex MUMPS. Figure 3. Open in new tabDownload slide Diagonal section through the FD impulse response in a homogeneous medium. (a) Splitting along the inline and crossline directions; (b) four-way splitting in the inline, crossline, and diagonal directions; (c) four-way splitting plus Li correction; (d) real MUMPS; (e) complex BICGSTAB; (f) complex MUMPS. 3.2. SEG/EAGE salt model Next, we use the SEG/EAGE salt model (Aminzadeh et al1995) to evaluate the performance of 3D migration without splitting in a complex inhomogeneous medium. The strong lateral variations between the sediments and the salt body and the high velocity of the salt are a challenge for the iterative methods. To filter out the spike velocity contrast used to simulate reflectors in the original model, we applied a median filter of size 7 × 7 × 7. The grid size, source position and source pulse are the same as for the homogeneous example. We computed FD impulses responses using BICGSTAB with a single term of the complex Padé expansion with rotation angles θ = 25° and θ = 45°. These impulse responses are compared to the FD impulse response computed using the direct solver MUMPS using three terms of the complex Padé response and θ = 25°. The higher rotation angles in the MUMPS solver tend to attenuate too much of the desired signal. Finally, we compare them to the 3D impulse response of reverse time migration using an FD scheme of the fourth order in space, second in time. For a single term of the complex Padé expansion, the MUMPS direct solver is two times slower than BICGSTAB. On the other hand, for higher-order complex Padé terms, the performance of BICGSTAB degrades rapidly, and increasing the angle θ does not help in this case. This agrees with the observation of Nichols (1997) that iterative methods have more difficulties in modelling the wavefield propagating at higher dip angles. Only the direct solver was able to compute the impulse response with three Padé terms. 3.2.1. Vertical sections Figures 4–6 show six vertical sections of the impulse response along the crossline direction (y) at inline positions x = 1500 m, x = 3880 m, and x = 4160 m, computed with 3D FD migration using BICGSTAB with a single complex Padé term with θ = 25° and θ = 45° and compare them to corresponding results obtained with the MUMPS solver with θ = 25°, two- and four-way splitting plus additional Li corrections with ten reference velocities at every six depth steps, and RTM. When comparing the parts of figures 4–6, four important observations can be made: (a) the images differ significantly in the shallow region. This is due to different artefacts, since because of its intrinsic restrictions, single-term FD migration cannot properly image events dipping above 45°. We see that the artefacts of the splitting+correction approach are much stronger than those of the iterative solver. The strongest differences are in the RTM image, since this method is, by design, the only one to treat these waves correctly. (b) The positioning of the migrated lower dipping events is very similar in all parts of these figures. The apparent changes in polarity are due to a slightly different positioning of the event in the inline (x) direction. (c) Apart from the different artefacts in the upper part, the Li correction causes a discontinuous behaviour of the wavefront. (d) The events of lower amplitude ahead of the main wavefront are not reproduced by the splitting+correction approach. Figure 4. Open in new tabDownload slide Sections of 3D FD migration impulse response in the plane x = 1500 m. From top to bottom: BICGSTAB, θ = 25°; BICGSTAB, θ = 45°; MUMPS, θ = 25°; two-way splitting plus Li correction; four-way splitting plus Li correction; RTM. Figure 4. Open in new tabDownload slide Sections of 3D FD migration impulse response in the plane x = 1500 m. From top to bottom: BICGSTAB, θ = 25°; BICGSTAB, θ = 45°; MUMPS, θ = 25°; two-way splitting plus Li correction; four-way splitting plus Li correction; RTM. Figure 5. Open in new tabDownload slide Sections of 3D FD migration impulse response in the plane x = 3380 m. From top to bottom: BICGSTAB, θ = 25°; BICGSTAB, θ = 45°; MUMPS, θ = 25°; two-way splitting plus Li correction; four-way splitting plus Li correction; RTM. Figure 5. Open in new tabDownload slide Sections of 3D FD migration impulse response in the plane x = 3380 m. From top to bottom: BICGSTAB, θ = 25°; BICGSTAB, θ = 45°; MUMPS, θ = 25°; two-way splitting plus Li correction; four-way splitting plus Li correction; RTM. Figure 6. Open in new tabDownload slide Sections of 3D FD migration impulse response in the plane x = 4160 m. From top to bottom: BICGSTAB, θ = 25°; BICGSTAB, θ = 45°; MUMPS, θ = 25°; two-way splitting plus Li correction; four-way splitting plus Li correction; RTM. Figure 6. Open in new tabDownload slide Sections of 3D FD migration impulse response in the plane x = 4160 m. From top to bottom: BICGSTAB, θ = 25°; BICGSTAB, θ = 45°; MUMPS, θ = 25°; two-way splitting plus Li correction; four-way splitting plus Li correction; RTM. Figures 7–9 show the corresponding comparison of vertical sections along the inline direction (x), at crossline positions y = 1500 m, y = 3880 m, and y = 4160 m. Again, note the good match among corresponding events dipping up to around 45°, but the differences for higher propagation angles, as well as regarding amplitudes and migration artefacts. As with the crossline sections, the strongest artefacts can be observed in the section at y = 1500 m. From these examples, we conclude that the splitting+correction approach has a number of drawbacks that can be avoided by the iterative solver at the expense of a higher computational cost. Figure 7. Open in new tabDownload slide Sections of 3D FD migration impulse response in the plane y = 1500 m. From top to bottom: BICGSTAB, θ = 25°; BICGSTAB, θ = 45°; MUMPs, θ = 25°; two-way splitting plus Li correction; four-way splitting plus Li correction; RTM. Figure 7. Open in new tabDownload slide Sections of 3D FD migration impulse response in the plane y = 1500 m. From top to bottom: BICGSTAB, θ = 25°; BICGSTAB, θ = 45°; MUMPs, θ = 25°; two-way splitting plus Li correction; four-way splitting plus Li correction; RTM. Figure 8. Open in new tabDownload slide Sections of 3D FD migration impulse response in the plane y = 3380 m. From top to bottom: BICGSTAB, θ = 25°; BICGSTAB, θ = 45°; MUMPS, θ = 25°; two-way splitting plus Li correction; four-way splitting plus Li correction; RTM. Figure 8. Open in new tabDownload slide Sections of 3D FD migration impulse response in the plane y = 3380 m. From top to bottom: BICGSTAB, θ = 25°; BICGSTAB, θ = 45°; MUMPS, θ = 25°; two-way splitting plus Li correction; four-way splitting plus Li correction; RTM. Figure 9. Open in new tabDownload slide Sections of 3D FD migration impulse response in the plane y = 4160 m. From top to bottom: BICGSTAB, θ = 25°; BICGSTAB, θ = 45°; MUMPS, θ = 25°; two-way splitting plus Li correction; four-way splitting plus Li correction; RTM. Figure 9. Open in new tabDownload slide Sections of 3D FD migration impulse response in the plane y = 4160 m. From top to bottom: BICGSTAB, θ = 25°; BICGSTAB, θ = 45°; MUMPS, θ = 25°; two-way splitting plus Li correction; four-way splitting plus Li correction; RTM. 3.2.2. Horizontal sections In the next set of figures, we compare horizontal sections of the BICGSTAB results to those of the other methods. As before, figures 10–13 show horizontal sections of the 3D impulse response computed by using six different methods: FD migration solved using BICGSTAB for one complex Padé term with θ = 25° (top left) and θ = 45° (top right), FD migration using the MUMPS direct solver with θ = 25° (centre left), two-way splitting plus Li correction (centre right), four-way splitting plus Li correction (bottom left), and RTM (bottom right). The section in figure 10 is at a very shallow depth level, where only waves with very high propagation angles should be present. This figure mainly illustrates the differences among the one-way methods regarding the migration artefacts at high propagation angles. Note that none of the one-way methods can be expected to correctly position these events. The only correct image at this shallow depth is obtained using RTM. Nonetheless, the main wavefront in the three-term MUMPS and splitting+correction solutions is slightly better positioned because of the 45° dip limit of the single-term complex Padé expansion used in BICGSTAB. Figure 10. Open in new tabDownload slide Sections of 3D FD migration impulse response in the plane z = 550 m. Top left: BICGSTAB, θ = 25°. Top right: BICGSTAB θ = 45°. Centre left: MUMPS, θ = 25°. Centre right: two-way splitting plus Li correction. Bottom left: four-way splitting plus Li correction. Bottom right: RTM. Figure 10. Open in new tabDownload slide Sections of 3D FD migration impulse response in the plane z = 550 m. Top left: BICGSTAB, θ = 25°. Top right: BICGSTAB θ = 45°. Centre left: MUMPS, θ = 25°. Centre right: two-way splitting plus Li correction. Bottom left: four-way splitting plus Li correction. Bottom right: RTM. Figure 11. Open in new tabDownload slide Sections of 3D FD migration impulse response in the plane z = 1050 m. Top left: BICGSTAB, θ = 25°. Top right: BICGSTAB θ = 45°. Centre left: MUMPS, θ = 25°. Centre right: two-way splitting plus Li correction. Bottom left: four-way splitting plus Li correction. Bottom right: RTM. Figure 11. Open in new tabDownload slide Sections of 3D FD migration impulse response in the plane z = 1050 m. Top left: BICGSTAB, θ = 25°. Top right: BICGSTAB θ = 45°. Centre left: MUMPS, θ = 25°. Centre right: two-way splitting plus Li correction. Bottom left: four-way splitting plus Li correction. Bottom right: RTM. Figure 12. Open in new tabDownload slide Sections of 3D FD migration impulse response in the plane z = 1350 m. Top left: BICGSTAB, θ = 25°. Top right: BICGSTAB θ = 45°. Centre left: MUMPS, θ = 25°. Centre right: two-way splitting plus Li correction. Bottom left: four-way splitting plus Li correction. Bottom right: RTM. Figure 12. Open in new tabDownload slide Sections of 3D FD migration impulse response in the plane z = 1350 m. Top left: BICGSTAB, θ = 25°. Top right: BICGSTAB θ = 45°. Centre left: MUMPS, θ = 25°. Centre right: two-way splitting plus Li correction. Bottom left: four-way splitting plus Li correction. Bottom right: RTM. Figure 13. Open in new tabDownload slide Sections of 3D FD migration impulse response in the plane z = 1550 m. Top left: BICGSTAB, θ = 25°. Top right: BICGSTAB θ = 45°. Centre left: MUMPS, θ = 25°. Centre right: two-way splitting plus Li correction. Bottom left: four-way splitting plus Li correction. Bottom right: RTM. Figure 13. Open in new tabDownload slide Sections of 3D FD migration impulse response in the plane z = 1550 m. Top left: BICGSTAB, θ = 25°. Top right: BICGSTAB θ = 45°. Centre left: MUMPS, θ = 25°. Centre right: two-way splitting plus Li correction. Bottom left: four-way splitting plus Li correction. Bottom right: RTM. Figures 11–13 demonstrate the quality of the one-term BICGSTAB solution. The impulse response is much closer to the theoretically better three-term MUMPS result than the one obtained with two-way splitting plus Li correction. Note again the noncausal artefacts of the latter approach ahead of the wavefront. From these figures, we see that the difference among corresponding events between MUMPS and BICGSTAB diminishes with depth, indicating that the iterative solver is sufficiently accurate for propagation directions below 45°. The differences between the 25° and 45° rotation angle BICGSTAB solutions are rather subtle. The main difference lies in the stronger attenuation of high-angle propagation for the 45° rotation angle. 3.2.3. Diagonal sections The final two figures 14 and 15 show the corresponding sets of images in the SW–NE and NW–SE diagonal directions. The most visible differences between the solutions is the high-frequency scatter in the real two and four-way splitting results. The results obtained using the complex Padé approximation contain far fewer migration artefacts. As before, the one-way wave-equation solutions differ much less between each other than from the RTM result in the last part. Figure 14. Open in new tabDownload slide SW-NE diagonal sections of 3D FD migration impulse response. From top to bottom: BICGSTAB, θ = 25°; BICGSTAB, θ = 45°; MUMPS, θ = 25°; two-way splitting plus Li correction; four-way splitting plus Li correction; RTM. Figure 14. Open in new tabDownload slide SW-NE diagonal sections of 3D FD migration impulse response. From top to bottom: BICGSTAB, θ = 25°; BICGSTAB, θ = 45°; MUMPS, θ = 25°; two-way splitting plus Li correction; four-way splitting plus Li correction; RTM. Figure 15. Open in new tabDownload slide NW-SE diagonal sections of 3D FD migration impulse response. From top to bottom: BICGSTAB, θ = 25°; BICGSTAB, θ = 45°; MUMPS, θ = 25°; two-way splitting plus Li correction; four-way splitting plus Li correction; RTM. Figure 15. Open in new tabDownload slide NW-SE diagonal sections of 3D FD migration impulse response. From top to bottom: BICGSTAB, θ = 25°; BICGSTAB, θ = 45°; MUMPS, θ = 25°; two-way splitting plus Li correction; four-way splitting plus Li correction; RTM. 4. Computation time When comparing splitting and nonsplitting solvers of 3D FD wave-equation migration, a word on computation time is indispensable. In our implementation, the BICGSTAB solution for one-term complex Padé FD migration took half as much time as the corresponding MUMPS solution and about five times as much as a corresponding solution based on two-way splitting plus Li correction with ten reference velocities. The three-term MUMPS solution consumed three times more time than the one-term solution. The three-term BICGSTAB solution converges in a much longer time than the three terms with the MUMPS solution. In other words, of the methods used in the above figures, BICGSTAB with a single term is intermediate between the six times more time-consuming three-term MUMPS solution and the five times less time-consuming splitting+correction approach. Our RTM implementation needed little less time than the MUMPS solution. Our results indicate that the combination of the complex Padé approximation with additional preconditioning may further reduce the computational cost, possibly bringing nonsplitting solutions within the range of more sophisticated corrections to splitting. Moreover, since low frequency is the main problem affecting the performance of BICGSTAB, multigrid techniques can certainly help to speed up these methods. Another possibility could be the combination of splitting and nonsplitting techniques for different frequency ranges or a hybrid approach combining iterative and direct solvers like the one suggested by Kao et al (1993). 5. Conclusions The performance of iterative methods for 3D downward continuation without splitting along inline and crossline directions improves markedly when using the complex Padé approximation. Numerical experiments using the SEG/EAGE salt model show the feasibility of FD migration without splitting for such a complex velocity model, being less than an order of magnitude slower than conventional splitting+correction approaches. We provided a simple analysis indicating that the complex Padé expansion improves the conditioning of FD downward continuation. For a single term in the Padé expansion, BICGSTAB outperforms the massively parallel direct solver by a factor of two. The performance of iterative methods degrades rapidly when more Padé terms are used but does not affect MUMPS in the same way. The computation time of MUMPS increases only linearly with the number of terms. The efficiency of full solutions without operator splitting still remains an issue. Our results indicate that the combination of the complex Padé approximation with additional preconditioning may further reduce the computational cost, possibly bringing nonsplitting solutions within the range of more sophisticated corrections to splitting. Moreover, since low frequency is the main problem affecting the performance of BICGSTAB, multigrid techniques can certainly help to speed up these methods. Another possibility could be the combination of splitting and nonsplitting techniques for different frequency ranges or a hybrid approach combining iterative and direct solvers. Acknowledgments We acknowledge the financial support from PETROBRAS and the sponsors of the Wave Inversion Technology (WIT) consortium, as well as the Brazilian Agencies Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES), Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) and Financiadora de Estudos e Projetos (FINEP). References Amazonas D , Costa J C , Schleicher J , Pestana R . , 2007 Wide-angle FD and FFD migration using complex Padé approximations , Geophysics , vol. 72 (pg. S215 - S220 ) 10.1190/1.2785813 Google Scholar Crossref Search ADS WorldCat Crossref Amestoy P R , Duff I S , Koster J , L'Excellent J Y . , 2001 A fully asynchronous multifrontal solver using distributed dynamic scheduling , SIAM J. Matrix Anal. Appl. , vol. 23 (pg. 15 - 41 ) 10.1137/S0895479899358194 Google Scholar Crossref Search ADS WorldCat Crossref Amestoy P R , Guermouche A , LExcellent J Y , Pralet S . , 2006 Hybrid scheduling for the parallel solution of linear systems , Parallel Comput. , vol. 32 (pg. 136 - 141 ) 10.1016/j.parco.2005.07.004 Google Scholar Crossref Search ADS WorldCat Crossref Aminzadeh F , Burkhard N , Kunz T , Nicoletis L , Rocca F . , 1995 3D modeling project: 3rd report , Learn. Edge , vol. 14 (pg. 125 - 128 ) 10.1190/1.1437102 Google Scholar Crossref Search ADS WorldCat Crossref Bamberger A , Engquist B , Halpern L , Joly P . , 1988 Higher order paraxial wave equation approximations in heterogeneous media , SIAM J. Appl. Math. , vol. 48 (pg. 129 - 154 ) 10.1137/0148006 Google Scholar Crossref Search ADS WorldCat Crossref Brown C L . , 1983 Application of operator separation in reflection seismology , Geophysics , vol. 48 (pg. 288 - 294 ) 10.1190/1.1441468 Google Scholar Crossref Search ADS WorldCat Crossref Claerbout J F . , 1985 , maging the Earth's Interior Oxford Blackwell Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Cole S . , 1989 Iterative methods for 3D finite-difference migration and modeling , Stanford Exploration Project Report SEP-60 pp 149–64 Costa J C , Mondini D , Schleicher J , Novais A . , 2011 A comparison of splitting techniques for 3D complex Padé Fourier finite difference migration , Int. J. Geophys. , vol. 2011 714781 10.1155/2011/714781 OpenURL Placeholder Text WorldCat Crossref Fei T , Etgen J . , 2002 Domain decomposition for 3D finite-difference depth extrapolation 72nd Annu. Int. Meeting SEG Technical Program Expanded Abstracts (pg. 1160 - 1163 ) 10.1190/1.1816855 Iserles A . , 1996 , First Course in the Numerical Analysis of Differential Equations Cambridge Cambridge University Press Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Kao J C , Li G , Yang C W . , 1993 Preconditioned iterative 3-D finite-difference depth migration or modeling on the CRAY T3D massively parallel processors 1993 SEG Annu. Meeting: SEG Technical Program Expanded Abstracts , vol. 12 (pg. 185 - 188 ) http://link.aip.org/link/?SEGEAB/12/185/1 10.1190/SEGEAB.12 Li Z . , 1991 Compensating finite-difference migration by multiway splitting , Geophysics , vol. 56 (pg. 1650 - 1660 ) 10.1190/1.1442975 Google Scholar Crossref Search ADS WorldCat Crossref Millinazzo F A , Zala C A , Brooke G H . , 1997 Square-root approximations for parabolic equation algorithms , J. Acoust. Soc. Am. , vol. 101 (pg. 760 - 766 ) 10.1121/1.418038 Google Scholar Crossref Search ADS WorldCat Crossref Mondini D , Costa J C , Schleicher J , Novais A . , 2012 Three-dimensional complex Padé FD migration: splitting and corrections , Int. J. Geophys. , vol. 2012 479492 10.1155/2012/479492 OpenURL Placeholder Text WorldCat Crossref Nichols D . , 1997 3D depth migration by a predictor-corrector method , Stanford Exploration Project Report SEP-70 pp 1–43 Rickett J , Claerbout J , Fomel S . , 1998 Implicit 3D depth migration by wavefield extrapolation with helical boundary conditions 68th Annu. Int. Meeting: SEG Technical Program Expanded Abstracts (pg. 1124 - 1127 ) 10.1190/1.1820086 Ristow D , Rühl T D . , 1997 3D implicit finite-difference migration by multiway splitting , Geophysics , vol. 62 (pg. 554 - 567 ) 10.1190/1.1444165 Google Scholar Crossref Search ADS WorldCat Crossref van der Vorst H A . , 1992 Bi-CGSTAB: a fast and smoothly converging variant of BI-CG for nonsymmetric linear systems , SIAM J. Sci. Stat. Comput. , vol. 13 (pg. 631 - 644 ) 10.1137/0913035 Google Scholar Crossref Search ADS WorldCat Crossref Wang Y . , 2001 ADI plus interpolation: accurate finite-difference solution to 3D paraxial wave equation , Geophys. Prospect. , vol. 49 (pg. 547 - 556 ) 10.1046/j.1365-2478.2001.00278.x Google Scholar Crossref Search ADS WorldCat Crossref Zhang G , Zhang Y , Zhou H . , 2000 Helical finite-difference schemes for 3D depth migration 70th Annu. Int. Meeting: SEG Technical Program Expanded Abstracts (pg. 862 - 865 ) 10.1190/1.1816209 Zhang L , Hua B , Calandra H . , 2005 3D Fourier finite difference anisotropic depth migration , SEG Technical Program Expanded Abstracts , vol. vol 24 (pg. 1914 - 1917 ) http://link.aip.org/link/?SEGEAB/24/1914/1 10.1190/SEGEAB.24 Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Crossref Zhang L , Rector J W , Hoversten G M . , 2003 Split-step complex Padé migration , J. Seismic Explor. , vol. 12 (pg. 229 - 236 ) OpenURL Placeholder Text WorldCat Zhang L , Rector J W , Hoversten G M , Fomel S . , 2004 Split-step complex Padé-Fourier depth migration , SEG Technical Program Expanded Abstracts , vol. vol 23 (pg. 989 - 992 ) http://link.aip.org/link/?SEGEAB/23/989/1 Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Zhang L , Rector J W , Hoversten G M , Fomel S . , 2007 Split-step complex Padé-Fourier depth migration , Geophys. J. Int. , vol. 171 (pg. 1308 - 1313 ) 10.1111/j.1365-246X.2007.03610.x Google Scholar Crossref Search ADS WorldCat Crossref © 2013 Sinopec Geophysical Research Institute
TI - Iterative methods for 3D implicit finite-difference migration using the complex Padé approximation
JF - Journal of Geophysics and Engineering
DO - 10.1088/1742-2132/10/4/045011
DA - 2013-08-01
UR - https://www.deepdyve.com/lp/oxford-university-press/iterative-methods-for-3d-implicit-finite-difference-migration-using-WZKz4zLvSj
VL - 10
IS - 4
DP - DeepDyve
ER -