Optimized Pseudo-Padé Fourier Migrator in Terms of Propagation Angles

A wide-angle Fourier migrator is proposed to more accurately image complex media with strong lateral velocity contrasts. The broadband wave-equation migrator is developed based on the pseudo-Padé approximation, where the Padé coefficients are independent of spatial coordinates, leading to a pure Fourier transform-based matching solution for one-way wavefield extrapolation. We use genetic algorithms to estimate the constant Padé coefficients more accurately than is feasible with conventional least-squares methods. Because of the global features of pure Fourier migrators, we present an angle-partitioning optimization scheme with dip focusing to improve the performance of the Fourier migrator for super-wide-angle waves and strong velocity contrasts. The wavefield gradient is used to calculate propagation angles during dual-domain wavefield extrapolation. Particular attention is paid to the first-order optimized pseudo-Padé Fourier (OPF1) migrator, which significantly improves the split-step Fourier (SSF) method for strong lateral variations at the cost of one additional Fourier transform in each step. Wavefield extrapolation based on the OPF1 method actually constitutes linear interpolation in the wavenumber domain between two split-step terms. We benchmark the OPF1 migrator with other typical migrators based on the exact dispersion equation. Numerical experiments with impulse responses, the SEG/EAGE salt model and 3D field data demonstrate the excellent performance and efficiency of seismic imaging with the OPF1 migrator.


I. INTRODUCTION
Fourier wave-equation migration is a rapidly developing area of research because of its computational efficiency in commercial applications. This approach has been widely recognized as a practical alternative to current industry techniques for accurate seismic imaging. A major impetus for studies on Fourier migrators is attributed to their many desirable properties, such as the analytical extrapolation of wavefields, preservation of amplitudes by honouring Snell's law, the simplicity of the algorithm and its high efficiency with fast Fourier transforms (FFTs), and its immunity to both grid dispersion and operator splitting errors.
The associate editor coordinating the review of this manuscript and approving it for publication was Nilanjan Dey.
Particularly, fast wavefield extrapolation with dual-domain implementation can be obtained by using a faster FFT (possibly with FFT chips). One of the key issues associated with Fourier migrators in handling wide-angle waves and lateral velocity variations simultaneously is the global features of the Fourier transform; this issue has led to the acknowledgement that velocity variations, propagation angles, and imaging accuracies are closely related at a variety of scales. However, high-order approximations involve extensive numerical calculations that do not conspicuously simplify the problem to be solved relative to other methods. Accordingly, in this article, we pursue a simple approximation and achieve an optimal trade-off among velocity variations, wide-angle waves, imaging accuracy, and computational efficiency. The motivation of this research is to obtain the highest imaging accuracy for wide-angle waves with the fewest FFTs.
Fourier imaging methods with different accuracies have been extensively studied in recent decades, beginning with the phase-shift propagator [1] for laterally homogeneous media and the split-step Fourier (SSF) method [2] and phasescreen propagator [3] for either weak-contrast media or small angles. Progress in this field can be classified into three categories: the approximation of the square-root operator, the coefficient optimization of Fourier propagators, and the partitioning of propagation angles. Slight modifications to phase-screen propagators have produced pseudo-screen propagators [3], [4], generalized-screen propagators (GSPs) [5], [6], and broadband constant-coefficient propagators [7] for moderate-contrast media. For high-contrast media with wide angles, the most appropriate approach may be the Rytov approximation [8] or separation-of-variables [9] screen propagators. These dual-domain propagators are generally formulated by the first-order approximation of the square-root operator, whereas higher-order approximations enhance the imaging accuracy but require additional Fourier transforms at each step, which considerably increases the computational time, especially for immense 3D cases. For high-contrast media with wide angles, currently, the most economical method might be the hybrid Fourier finitedifference (FFD) method [10], which is not a pure Fourier migrator and hence needs an additional implicit FD correction. It is worth mentioning that Zhang et al. [11] applied firsttype Chebyshev polynomials to economize the result of the Taylor expansion of the square-root operator. The resulting Chebyshev Fourier migrator can be optimized to handle highcontrast media with wide angles but at the cost of four Fourier transforms in each step (twice the number of transforms in the SSF method). In this study, we pursue a simple rational approximation to the square-root operator; the resulting pure Fourier migrator needs only three Fourier transforms in each step for high-contrast media with wide angles.
The second category of Fourier migrators results from the optimization of constant expansion coefficients [12], [13] by least-square methods. Most optimization schemes used for the FFD method [10] can be directly applied to Fourier migrators, but limited improvement is achieved because the aforementioned Fourier migrators are based mainly on the low-accuracy Taylor expansion. The third type of scheme used to improve the SSF method attempts to directly address the difficulty associated with the global features of Fourier migrators; such schemes include the phase shift plus interpolation (PSPI) [14], windowed phase-screen [15], and optimum separable approximation [16] methods. These techniques can image steep dips in complex media with strong velocity contrasts but require many more FFTs or additional computations. In this article, we take advantage of the genetic algorithm (GA) for the global optimization of Fourier expansion coefficients and propose an angle-partitioning optimization method to enhance the localization of Fourier migrators for individual propagation angles.
Most rational functions converge faster than the Taylor expansion because of their well-posed approximation behaviour for low-order terms [17], [18]. Most Fourier migrators obtained by the poorly convergent Taylor expansion are characterized by few benefits, even with many highorder terms. In contrast, rational approximations with superior accuracy can significantly improve the performance of Fourier migrators for constructing the square-root operator. In general, a first-order or, at most, a second-order approximation is adequate for common seismic imaging applications in high-contrast media with wide angles. For example, the Padé approximation, an efficient rational function, has been used to construct various hybrid Fourier matching solutions for one-way propagation problems in high-contrast media; the corresponding migrators include the Padé FD migrator [19], the split-step Padé propagator [20], and the FFD migrator [10]. These rational approximations to the square-root operator lead to variable Padé coefficients as functions of the corresponding spatial coordinates and therefore require an additional implicit FD implementation.
To avoid additional implicit FD implementations, in this article, we propose a pseudo-Padé expansion to the squareroot operator; this approach has been used for the Bloch-Horowitz effective Hamiltonian [21]. The resultant Padé coefficients are constant, allowing a pure Fourier transformbased matching solution for seismic imaging. In particular, the first-order pseudo-Padé Fourier migrator significantly improves the SSF method for strong lateral variations at the cost of one additional Fourier transform in each step. Wavefield extrapolation is actually a linear interpolation problem in the wavenumber domain between two split-step terms. Thus, determining how to optimize the constant Padé coefficients in terms of the maximum velocity contrasts and designed propagation angles is crucial to assure the fast convergence of a low-order pseudo-Padé expansion. We use the GA to optimize the constant Padé coefficients for all velocity contrasts. To localize the broadband pseudo-Padé Fourier migrator for individual angles, we apply an anglepartitioning optimization method [22] for dip focusing to improve the wide-angle accuracy of the pseudo-Padé Fourier migrator.
In this paper, a novel pure Fourier wave-equation migrator is proposed. The rest of this paper is organized as follows. In Section II, the proposed migrator is clearly described, and we demonstrate its effectiveness by performing relative phase-error analyses. In Section III, we conduct some numerical experiments to illustrate the proposed scheme. Finally, in Section IV, we discuss the practicality of the proposed migrator and summarize its main advantages.

A. FOURIER MIGRATOR USING THE PSEUDO-PADÉ APPROXIMATION
We slice heterogeneous media horizontally into a stack of slabs. Let u(k x , z) represent the 2D time-harmonic scalar VOLUME 8, 2020 wavefield in the frequency-wavenumber domain, where z is the depth and k x is the wavenumber with respect to the x-coordinate. For a laterally heterogeneous slab spanning from z to z + z, Fourier wavefield extrapolation can be generally expressed in the frequency-wavenumber domain as where k 2 x + k 2 z = k 2 o and the reference wavenumber is k 0 = ω/υ 0 with the reference velocity υ 0 for the slab. The intermediate wavefieldû(k x , z) is given by different Fourier migrators for different accuracies.
As described in Fu [23], applying the plane-wave representation of the Hankel function to the generalized Lippmann-Schwinger integral equation [24] leads to the following approximation of the square-root operator (refer to Appendix A for details): where n is the refractive index and the normalized wavenumbers arek x = k x k 0 andk z = k z k 0 . Because the Born approximation is used to reduce the two-way Lippmann-Schwinger integral equation to a one-way version, (2) can be regarded as the Born dispersion equation with an accuracy curve that is much better than that of the SSF method [23]. However, wavefield extrapolation by (2) using Fourier transforms is unstable because whenk x ≈ 1 (90 • propagation angle), the denominator 1 −k 2 x is close to 0. To avoid a singularity, we can apply the Padé approximation to (2), resulting in stable and well-posed wavefield extrapolation. However, the Padé operator is 'local' in that the Padé coefficients vary with both n andk x . Hence, the implementation of the differential operator requires the FD scheme. We rewrite (2) as To fit (3) into the pure Fourier matching framework, we propose the following pseudo-Padé expansion (refer to Appendix B for further details): where the pseudo-Padé coefficients a j and b j are independent of the refractive index n. We substitute (4) into (3) to yield the pseudo-Padé dispersion equation: This equation theoretically has no singularity and can be implemented for unconditionally stable wavefield extrapolation. Formulating the intermediate wavefieldû(k x , z) by (5) yields the following Fourier wavefield extrapolation (1), where we set C j = a jk Since the well-posed rational approximation in (5) is characteristic of the fast convergence of low-order terms, the firstorder equation is usually adequate for common seismic imaging applications. Thus, we take the first-order approximation, and (6) becomes We see that wavefield extrapolation by (7) actually constitutes linear interpolation in the wavenumber domain between two split-step terms. Therefore, the first-order pseudo-Padé Fourier migrator significantly improves the SSF method for strong lateral variations at the cost of only one additional Fourier transform in each step. Note that wavefield extrapolation by (7) tends to be numerically unstable unless some restrictions are imposed on C j . That is, the wavefield energy for some n values will increase during extrapolation. From (6), a stable numerical implementation requires This criterion is not convenient because k 0 , n, and z are involved. It follows from (8) that This inequality should be used as a constraint in the optimization procedure employed to search for the optimal a j and b j . Wavefield extrapolation by (6) under condition (10) will be unconditionally stable.

B. OPTIMIZATION OF COEFFICIENTS BY GENETIC ALGORITHMS
The optimization procedure employed to search for the optimal a j and b j can be defined by minimizing the following cost function: where θ is the propagation angle, which varies from 0 • to the designed maximum angle ϕ. From (5), the dispersion error of the pseudo-Padé approximation can be written as We see that the optimization procedure based on (11) and (10) is a constrained nonlinear optimization problem for multiparameter optimization. In general, local optimization techniques such as the least-squares method and its variants rely on exploiting limited information derived from a comparatively small number of models. Thus, a local optimization method is likely to be trapped in local minima if the starting values (a j and b j ) are too far from the true values [25], [26]. Theoretically, local optimization methods have difficulty searching for the best fit over a wide range of angles.
As an alternative, GAs are nonlinear global optimization techniques based on the natural processes of biological evolution and have been applied extensively to geophysical optimization problems, e.g., seismic inversion [27], attribute selection [28], migration velocity estimation [29], and residual statics corrections [30]. GAs are driven completely by stochastic means. Essentially, the search mechanism does not follow a deterministic set of rules, giving GAs the ability to escape local minima without the need for starting values. Thus, we use GAs to handle the optimization problem and search for the optimal a j and b j over a wide range of angles.
The pseudo-Padé Fourier migrator proposed in this article is a constant-coefficient operator. Because of its global features, no single values of a j and b j can adapt to all the velocity contrasts over a wide range of angles. To increase the localization of the constant-coefficient operator, we use angle-partitioning optimization [22] for dip focusing. The coefficients a j and b j are estimated in the GA optimization procedure for both the small-angle range (0 • ∼50 • ) and the wide-angle range (50 • ∼75 • ). The resultant coefficients for the first-order optimized pseudo-Padé Fourier (OPF1) migrator are a 1 = −0.627 and b 1 = 0.122 over 0 • ∼50 • and a 1 = −0.654 and b 1 = 0.087 over 50 • ∼75 • , with the phase-error analyses at refractive index n = 0.7 depicted in Figure 1 at different relative phase errors. We see that the small-angle curve (solid line) varies within 0 • ∼50 • under a relative phase error of 1% but allows a maximum propagation angle reaching 58 • under a relative phase error of 5%. In contrast, the wide-angle curve (dashed line) varies within 0 • ∼57 • under a relative phase error of 2% but allows a maximum propagation angle reaching 62 • under a relative phase error of 5%. An analysis of Figure 1 reveals that angle-partitioning optimization has the ability to improve the performance of the OPF1 migrator and further improve both the small-angle curve and the wide-angle curve at individual propagation angles for a particular target area. It is worth noting that the OPF1 migrator for strong velocity contrasts still preserves a dual-domain algorithmic structure that uses three Fourier transforms alone to advance the wavefield.

C. OPTIMIZATION BY ANGLE PARTITIONING
In equation (5), n andk x are cross-coupled. Thus, the OPF1 migrator is implemented in the frequencywavenumber domain. To perform the angle-partitioning optimization procedure with the OPF1 migrator, we need to calculate the wave propagation angle at any given point in space in the frequency-wavenumber domain; this approach has been widely used in seismic migration and velocity analysis [31]- [33]. In general, there are two kinds of calculation schemes: ray-based [34], [35] and wavefield gradient methods [22]. Based on the corresponding high-frequency limitation, ray-based methods are easier to implement and have a much lower computational cost. However, the results of raybased methods are difficult to implement with Fourier-based wavefield extrapolation methods. In contrast, the wavefield gradient can be easily obtained from a Fourier-based one-way propagator, and the wave propagation angles can also be calculated at any specified point in space for a single-frequency component. In addition, the wavefield gradient method is fast and easy to implement with a Fourier wave-equation migrator. Thus, we adapt the wavefield gradient method to handle the angle-partitioning optimization procedure.
Based on the wavefield gradient vector, we obtain the wave propagation angle with respect to the z direction: where ∂u(x,z) ∂x and ∂u(x,z) ∂z are the wavefield gradients with respect to the x-and z-coordinates, respectively. Then, we use the OPF1 migrator for the calculation of the wave VOLUME 8, 2020 propagation angle. Equation (7) can be expressed as (14) where FT x and FT −1 x denote the forward and inverse Fourier transforms, respectively, from the space domain to the wavenumber domain in the x direction. Based on (14), the wavefield gradient along the x-coordinate can be written as Based on the unintegrated form of (14), the wavefield gradient in the z direction is expressed as Because the OPF1 migrator is a dual-domain propagator, the wavefield gradient at a given point in space is usually a complex number. Here, we choose the absolute values of (13) as the wave propagation angles. Based on (13), (15), and (16), we have the ability to calculate the wave propagation angle at any given point in the frequency-wavenumber domain.   3) Use a GA to search for the optimal coefficients a 1 and b 1 for the small-angle range and the wide-angle range, respectively. 4) Store these constant pseudo-Padé coefficients in a lookup table. 5) Assign the stored coefficients a 1 and b 1 for each frequency component at every point in space and perform wavefield extrapolation with (7).

D. RELATIVE PHASE-ERROR ANALYSES
In this section, we will discuss the accuracy of the OPF1 migrator with applications involving angular spectra and dispersion circles. The optimized Chebyshev Fourier (OCF) method [11], which is a pure Fourier migrator, has sufficiently high accuracy (approximately 60 • at a relative error of 1%) for almost all velocity contrasts. Thus, we select the OCF method as a reference to assess the OPF1 migrator. Figure 3 shows the angular spectra of various methods for a relative phase error of 1%. A large refractive index (n) denotes weak-contrast media, while a small value denotes high-contrast media. Obviously, the SSF method displays the worst performance among all the methods tested. The first to fourth orders of the GSPs exhibit different improvements for weak contrasts but rapidly degenerate to the performance of the SSF method with a decrease in the refractive index. We see that the FFD, OPF1, and OCF methods perform well for almost all the lateral velocity variations, and each curve is nearly horizontal. The OPF1 method displays a significant improvement over the GSP and SSF methods and almost approaches the OCF method in terms of accuracy. Although we can see that the OCF method yields a slight improvement over the OPF1 method, the OCF migrator needs one additional Fourier transform in each time step. To assess the global properties of the OPF1 migrator, we evaluate the proposed method by comparing its dispersion circles with those of the SSF, FFD, and exact dispersion relations for refractive index values of n = 0.3 (strong contrast), n = 0.65 (moderate contrast), and n = 0.9 (weak contrast) (Figure 4). Among these methods, the SSF method has the lowest accuracy for all the n values tested. For n = 0.9, the dispersion circles of the OPF1 and FFD methods fit the exact dispersion relations at almost all angles of wave propagation. For n = 0.65 and n = 0.3, the dispersion circles of the OPF1 and FFD methods have high accuracy at small angles, and the accuracy decreases with increasing angle. Obviously, the maximum propagation angles of the OPF1 method exceed those of the FFD method for almost all velocity contrasts, except in areas with very weak lateral velocity contrasts (Figure 3 and Figure 4).

III. EXAMPLES
In this section, we use impulse responses, the 2D Society of Exploration Geophysicists (SEG)/European Association of Geoscientists and Engineers (EAGE) salt model, and real 3D field data to illustrate the advantages of the OPF1 method. All of the examples are implemented in the Seismic Unix (SU) software system and are simulated on a rack-mounted server with a 12-core Intel R Xeon R Gold 5118.

A. IMPULSE RESPONSE
To illustrate the performance of the OPF1 method, we first calculate several impulse responses. A 2D homogeneous medium is defined for a system of 256×128 grids with a grid spacing of 10 m. The real velocity is V = 4000 m/s. A Ricker wavelet with a dominant frequency of 35 Hz is used as the source and is located at the centre of the upper model boundary. The travel time is 250 ms, and the sampling interval is 2 ms. The frequency range is from 1 Hz to 100 Hz with 101 components. Figure 5 shows vertical slices of the snapshots at 250 ms. These plots show the impulse responses obtained by the SSF and OPF1 methods. The refractive indexes used for each plot in Figure 5 are (a) n = 0.35 (strong contrast), (b) n = 0.65 (moderate contrast), and (c) n = 0.9 (weak contrast). The dashed semicircles are the exact impulse responses and serve as the quality benchmark. Obviously, the impulse responses of the SSF method deviate significantly from the theoretical curves except in the weak-contrast media (n = 0.9). In contrast, the OPF1 method always shows a stable performance, and the results remain close to the theoretical positions within a propagation angle of 60 • for all three n values tested. Figure 6 shows the superpositions of the vertical slices of the impulse responses obtained by the FFD (left part of each plot) and OPF1 (right part of each plot) methods. When using a grid interval of 10 m, the OPF1 method shows VOLUME 8, 2020 a slight improvement over the FFD method for n = 0.65 (Figure 6a). When a relatively coarse grid of 20 m is used, the FFD method produces significant numerical dispersion and strongly distorts the impulse responses at wide angles, as shown in Figure 6b. In contrast, the OPF1 method, which is immune to numerical dispersion, always produces a good result with high accuracy and no numerical dispersion.

B. 2D SEG/EAGE SALT MODEL
In the second example, we compare the 2D migrated images obtained using the OPF1 migrator with those acquired using the SSF, FFD and PSPI methods. These methods have different imaging accuracies and computational efficiencies; therefore, they constitute good benchmarks for testing the novel OPF1 migrator. As mentioned above, we conclude that the main advantage of the OPF1 migrator is that it can be applied to image steeply dipping reflectors with strong velocity contrasts. The SEG/EAGE salt model (Figure 7a) has steep dip reflectors (the largest value is approximately 70 • ) and extreme lateral velocity perturbations (n min ≈ 0.35). Thus, we select the SEG/EAGE salt model to demonstrate the validity of the OPF1 method. Figure 7b shows the calculation results for sinθ (θ is the wave propagation angle with respect to the z direction at a given point) at 60 Hz based on (13), (15) and (16). Obviously, waves travel though the salt body at wide propagation angles, and the salt body complicates the subsalt wavefield. According to the characteristics of the salt model, we set the boundary angle to be θ b = 50 • . For θ (x, z) < 50 • , the given point (x, z) falls within the small-angle range (the value of sinθ is 0 in Figure 7c, 7d), whereas it falls within the wide-angle range otherwise (the value is 1 in Figure 7c, 7d). We note that the wide-angle parts for 30 Hz and 35 Hz (Figure 7c, 7d) are very similar overall. In this case, we can calculate the wave propagation angles for every five frequency components and compute the others via interpolation; this approach is sufficient for ensuring the accuracy of the OPF1 migrator and can improve the computational efficiency. Once the wave propagation angles are obtained, we can assign the stored coefficients for all the frequency components involved at every point in space. To further improve the accuracy of the OPF1 migrator, we can calculate additional single-frequency components and perform fewer interpolations.
The salt model contains several key targets (illustrated in Figure 7a) to demonstrate the validity of seismic migration methods, such as a strong velocity contrast along steep salt flanks (A), horizontal subsalt interfaces (B) and steep subsalt faults (C and D). In Figure 8, we show the resulting poststack depth-migrated images for the salt model from the SSF, FFD, PSPI, and OPF1 methods. From the outline of the salt body (marked by the blue lines in Figure 8), we note that the salt body is well imaged by almost all the migration methods tested, except for the SSF method. The SSF method works accurately only for either weak contrasts or small propagation angles, and it leads to a poor image of the steep salt root (Figure 8a). The FFD method, which involves one implicit FD correction and two Fourier transforms for each step, is not a pure Fourier method; Figure 8b reveals that the primary structural elements, except the steep fault C, are imaged well by the FFD method with weak migration noise. The PSPI method, which involves multiple Fourier transforms for each step, is implemented using wavefield interpolation. In Figure 8c, we see that the PSPI method produces a good image with high accuracy and weak migration noise; similar to the result with the FFD method, however, the steep subsalt fault C is not well imaged. The steep subsalt fault C is difficult to image for the above migration methods because the highvelocity salt body acts as an acoustic lens, spreading energy in a 'random' way [36]. In addition, the PSPI method attenuates the boundary reflections of the steep salt flank A; this can be explained by the migrator being influenced by attenuation at high wavenumbers. The OPF1 method, which requires only three Fourier transforms for each step, is a pure Fourier method. Compared with the FFD method, the OPF1 migrator partly attenuates coherent noise in the subsalt region because it provides better focusing capabilities for wideangle waves. Moreover, compared with the PSPI method, the OPF1 method produces better images of the steep salt root and the complex subsalt structures, particularly for the steep subsalt fault C (marked by the red arrow in Figure 8d), which cannot be imaged by most of the other methods (Figure 8d). Thus, these migration results demonstrate the validity of the OPF1 method in imaging complex structures with strong velocity contrasts. VOLUME 8, 2020

C. FIELD DATA
To further test the imaging capability of the OPF1 method in commercial applications, we consider real 3D field data acquired in eastern China that contain a number of complex geological structures. The complexities encountered during seismic imaging in this complex region are summarized as follows: (1) small faults/blocks and thin-bedded depositional units lead to low signal-to-noise (S/N) ratios and strong lateral velocity perturbations, and (2) the data set contains several steep faults and dipping interfaces. Conventional seismic migration methods have difficulty describing wideangle waves with strong lateral velocity variations. Hence, in Figure 9, we compare the results of pre-stack depth migration with the real data set using commercial software (Kirchhoff) and the OPF1 method. We easily discover that the OPF1 migrator produces a higher-quality image of the dipping interfaces and small faults/blocks, which are imaged in place, than the conventional Kirchhoff depth migration approach for this real dataset (marked by the blue circles in Figure 9). In the depth-migrated section using the Kirchhoff method (Figure 9a), the small faults/blocks in deep areas cannot be characterized due to a poor imaging quality. In addition, steeply dipping events are not well reconstructed and are smeared by coherent noise from incorrectly migrated data. The OPF1 migrator (Figure 9b) not only greatly improves the imaging quality in both shallow and deep areas, leading to a sharp image comprising small faults/blocks and steep fault planes, but also focuses steeply dipping events and suppresses coherent noise. The improvements in this migration result (Figure 9b) are consistent with the observations obtained from the SEG/EAGE salt model test, as illustrated in Figure 8. Obviously, the OPF1 migrator is more proficient at focus the energy of high-wavenumber fields than conventional migration methods; thus, the proposed approach can be recognized as a practical alternative to current commercial migration techniques.

IV. DISCUSSION AND CONCLUSION
Both relative phase-error analyses and numerical experiments demonstrate that the OPF1 method can image steep dips for almost all velocity contrasts. The maximum propagation angles of the OCF migrator [11] are always approximately 60 • under a relative error of 1%, and these values are approximately 5 • higher than those of the OPF1 method ( Figure 3). However, four Fourier transforms are required for each depth slice to advance the wavefield (Table 1). In contrast, the coefficients of the OPF1 method, which are calculated before migration, are stored in a lookup table (Figure 2), and only three Fourier transforms are required for the OPF1 method (Table 1); thus, the computational efficiency of the OPF1 method is higher than that of the OCF method. In general, the maximum propagation angle of 55 • (for a relative phase error of 1%) is adequate for common seismic imaging applications. Compared with the SSF method, the OPF1 migrator can handle wide-angle waves and strong lateral velocity variations simultaneously at the cost of one additional Fourier transform for each depth slice (Table 1). Therefore, we believe that the OPF1 method achieves an acceptable trade-off between the migration accuracy and computational efficiency.
To handle wide-angle waves and strong lateral velocity contrasts simultaneously, we propose an optimized pseudo-Padé Fourier migrator using GAs and an angle-partitioning optimization scheme. In this paper, we first apply a pseudo-Padé expansion to the square-root operator, and the resultant coefficients are constant, leading to a pure Fourier transformbased matching solution for seismic imaging. Then, we use a GA to search for the optimal constant Padé coefficients. Finally, we present an angle-partitioning optimization scheme for dip focusing to reduce the phase error for wide angles in high-contrast media. In this case, the wavefield gradient is used to calculate propagation angles. Based on our relative phase-error analyses and numerical experiments, the main advantages of the OPF1 method can be summarized as follows: 1) The OPF1 migrator is adequate for common one-way wavefield extrapolation applications, and its accurate dip angles (approximately 55 • at a relative phase error of 1%) are similar to those of the FFD and OCF methods for arbitrary lateral velocity contrasts.
2) The OPF1 migrator requires only three Fourier transforms for migration, and it is obviously superior to most of the current Fourier methods in terms of the computational cost. In addition, this approach has the potential ability to improve the computational efficiency of the transform by using a much faster FFT (possibly with FFT chips) compared to the traditional FFD method.
3) The OPF1 migrator is a pure Fourier method; thus, it is immune to numerical dispersion for coarse grids and high-frequency waves. In addition, this migrator has no operator splitting errors for 3D cases; thus, it can be easily extended from two to three dimensions.

A. BORN DISPERSION RELATIONSHIP FOR THE ONE-WAY LIPPMANN-SCHWINGER INTEGRAL EQUATION
This appendix illustrates a brief derivation of (2). Wave propagation in heterogeneous media can be described by the generalized Lippmann-Schwinger integral equation. We slice heterogeneous media horizontally into a stack of heterogeneous slabs and apply the plane-wave representation of the Hankel function to the generalized Lippmann-Schwinger integral equation. We obtain the following one-way Lippmann-Schwinger integral equation: where u is the seismic displacement vector, k x and k z are the horizontal and vertical wavenumbers, respectively, k z is the vertical wavenumber related to the media immediately following the slab, k 0 is the reference wavenumber, F x (k x , z) = FT x {ik 0 [n (r) − 1] u (r)} with FT x being the Fourier transform from x to k x , and n (r) is the refractive index, with r being the position vector. This equation preserves the main characteristics of wavefields in heterogeneous media and reduces the corresponding computational cost.
We normalize the wavenumbersk x = k x k 0 ,k z = k z k 0 ,k x = k x k 0 , andk z = k z k 0 , where k 0 is the background wavenumber of the adjacent media. For convenience, we take the refractive index n(r) of the slab as different n values to avoid the convolution of F x (k x , z + z); thus, the following analysis within the Fourier domain proceeds by decomposing the heterogeneous slab into a series of homogeneous slabs. In this sense, (A-1) can be written as where O(n) is the relative slowness perturbation defined as O(n) = n(r) 2 − 1, with n(r) being the refractive index. (A-2) can be written as With the Taylor series expansion of the arctangent terms of (A-3) inside the bracket for and taking the first-order term, we obtain From equation A-6, we obtain the following degenerate dispersion equation for the one-way Lippmann-Schwinger integral equation for heterogeneous slabs: Letk x ≈k x without considering refracted waves, and substitute this relation into (A-7) to yield the following degenerate dispersion equation:

B. PSEUDO-PADÉ APPROXIMATION
This appendix illustrates a brief derivation of (4). The Fourier matching algorithms based on (3) may not be sufficiently accurate for forward wave propagation in complex media with strong velocity contrasts. To avoid additional VOLUME 8, 2020 implicit FD implementation steps, we consider the following separation-of-variables operator decomposition: In this approach, we need to address the following two problems: the construction of the splitting operators f j k x and g j (n) and the application of the Fourier transform algorithm to (B-1). Initially, from the physical understanding of (B-1), the best way to construct the splitting operators f j k x and g j (n) might be to take g j (n) as a function of the split-step term (i.e., g j (n) = (n − 1) j ) and then to approximate f j (k x ) as a term of the rational functions; that is, where the coefficients a j and b j do not vary with the refractive index n and can be modified to reduce dispersion errors. Actually, because |k x | ≤ 1 for one-way propagation, the term (1 −k 2 x ) −1 / 2 can be approximated by the following pseudo-Padé expansion: From 2014 to 2015, he was a Visiting Researcher with Rice University. He is currently a Professor with the School of Geosciences, China University of Petroleum (East China), Qingdao, China. His research interests include seismic data processing, leastsquares reverse time migration, Gaussian beam imaging, and full waveform inversion. VOLUME 8, 2020