An optimised one-way wave migration method for complex media with arbitrarily varying velocity based on eigen-decomposition

The classical one-way generalised screen propagator (GSP) and Fourier finite-difference (FFD) schemes have limitations in imaging large angles in complex media with substantial lateral variations in wave velocity. Some improvements to the classical one-way wave scheme have been proposed with optimised methods. However, the performance of these methods in imaging complex media remains unsatisfying. To overcome this issue, a new strategy for wavefield extrapolation based on the eigenvalue and eigenvector decomposition of the Helmholtz operator is presented herein. In this method, the square root operator is calculated after the decomposition of the Helmholtz operator at the product of the eigenvalues and eigenvectors. Then, Euler transformation is applied using the best polynomial approximation of the trigonometric function based on the infinite norm, and the propagator for one-way wave migration is calculated. According to this scheme, a one-way operator can be computed more accurately with a lower-order expansion. The imaging performance of this scheme was compared with that of the classical GSP, FFD and the recently developed full-wave-equation depth migration (FWDM) schemes. The impulse responses in media with arbitrary velocity inhomogeneity demonstrate that the proposed migration scheme performs better at large angles than the classical GSP scheme. The wavefronts calculated in the dipping and salt dome models illustrate that this scheme can provide a precise wavefield calculation. The applications of the Marmousi model further demonstrate that the proposed approach can achieve better-migrated results in imaging small-scale and complex structures, especially in media with steep-dipping faults.


Introduction
The pre-stack depth migration can better deal with lateral velocity changes compared with the pre-stack time migration. The one-way wave depth migration (OWDM) has played a vital role in seismic migration in past decades. Because the one-way wave equation is separated from the full-wave equation, the OWDM method inherits the ad-vantages of the correct phase information of the two-way wave method and compensates for the shortcomings of the ray migration method. Therefore, the OWDM method has been widely used for collecting practical seismic information. Extensive research has been carried out on oneway wave imaging algorithms. Approximation expressions have been implemented extensively to improve the OWDM method. The initial phase-shift method, proposed to solve the wave equation in the frequency-wavenumber domain, is limited by the geophysical assumption of a homogeneous media (Gazdag 1978). This assumption is too stringent as actual media feature significant lateral changes in velocity. To overcome this issue, some improved schemes have been derived, such as the phase-shift plus interpolation scheme and the subsequently developed split-step Fourier (SSF) scheme (Gazdag & Sguazzero 1984;Stoffa et al. 1990). Because of the limited compensation for low-order perturbation, these schemes can only offer correct imaging when the media present mild lateral velocity variations. To further improve the ability to image medium with significant lateral velocity changes and large-angle structures, the Fourier finite-difference (FFD) and almost simultaneously proposed generalised screen propagator (GSP) schemes have started serving seismic data imaging (Ristow & Rühl 1994;Wu 1994;Jin & Wu 1999;De Hoop et al. 2000;Huang & Fehler 2000;Le Rousseau & De Hoop 2001). Both methods use a relatively higher-order correction term to process the lateral velocity variability. Consequently, better imaging results can be obtained in imaging steep dip angles than the SSF method. Cheng et al. (2001) proposed the finite-difference depth migration method in the frequency-space domain to better adapt velocity variations. Various propagators with distinct precision and efficiency can be produced by different approximation strategies, such as pseudo-screen operators, complex screen operators, windowed screen propagators, optimised product propagators and optimal approximate propagators. These methods, based on the decomposition of the velocity into background and disturbance velocities for wavefield extrapolation, lead to better results compared to the SSF method, yet they remain unable to correctly image the wavefield in complex media.
Based on the global optimum theory, some scholars have attempted to use more accurate multinomial coefficients in approximate expansions to increase the imaging accuracy of complex media with large angles, and developed optimisation strategies for the OWDM method (Biondi 2002;Zhu et al. 2008;Jia & Wu 2009;Chen 2010;Sun et al. 2010;Zhang et al. 2010). By using methods for approximating the propagator of the one-way wave equation, these approaches produced several improvements in large-angle imaging. However, they did not eliminate the limitation of the classical multistep calculation method in a mixed domain to suit for arbitrarily variant velocity. Therefore, a new method of using approximation theory to approach the correct velocity in the one-way wave equation needs to be further developed.
Eigen-decomposition is a mathematical method generally used for seismic data processing and inversion (Luo et al. 2003;Luk et al. 2005;Gheimasi et al. 2010;Wen et al. 2017). Some studies used this method to perform wave extrapolation in OWDM, but the performance was still inadequate. Kosloff & Kessler (1987) proposed an accurate depth migration method using eigendecomposition to improve the conventional phase-shift method. Yan & He (1994) derived a seismic wave propagating matrix in two-dimensional transversely isotropic media using the eigen-decomposition method and pointed out that this scheme can yield an accurate and straight forward propagator. Grimbergen et al. (1998) proposed an eigen-decomposition extrapolation method based on modal expansion, which divided the wavefield into guided and radiated wavefields. You et al. (2018) used the eigendecomposition method to achieve the amplitude-preserved calculation in OWDM. You et al. (2019) developed a scheme using matrix multiplication for wavefield extrapolation in laterally varying media, and obtained similar calculation results of the Helmholtz operator with the eigen-decomposition scheme. However, these methods are limited by their finite calculation accuracy unless a higher-order approximation with extensive computation is used. Trigonometric function approximations based on the Legendre, Taylor series and Chebyshev polynomial expansions have been used in physics, seismic modeling, and data processing and inversion (Dikhaminjia et al. 2015;Naghadeh & Morley 2017;Araujo & Pestana 2019). Nevertheless, few studies exploited the infinite-norm approximation of the trigonometric function to improve seismic wave exploration and imaging.
In this study, a strategy for eigenvalue and eigenvector decomposition of the Helmholtz matrix was investigated for OWDM. A novel scheme was proposed for computing the one-way propagation operator with the best polynomial approximation based on the infinite norm. The computing accuracy of the Taylor series and the least-squares (LS) approximation were used to compare with our scheme. The imaging performance was verified in the impulse response of three typical models in the imaging of media with significant lateral inhomogeneity and complex structures. The migrated sections computed by the classical one-way GSP, staggered-grid finite-difference (SFD) methods were compared with those obtained by the proposed scheme. Meanwhile, full-wave-equation depth migration (FWDM) methods deal with the wavefield with less approximation for the wave equation and can imaging waves at large angles compared to the OWDM methods (Kosloff & Baysal 1983;Sandberg & Beylkin 2009;You et al. 2016). The method proposed in You et al. (2016) overcame the surface constant velocity assumption and obtained a better imaging amplitude than the conventional RTM method. In their study, the Helmholtz operator was used to achieve the wavefield exploration based on the first-order stress-velocity acoustic equation in matrix form, which relates to the proposed scheme. Therefore, we adopted this method as a reference to verify the proposed scheme. Finally, to further demonstrate the ability of the proposed scheme for imaging complex media with dipping angles, detailed comparisons with the Marmousi model were carried out. 777 2. Depth extrapolation with the best approximation algorithm

One-way wave-equation depth extrapolation
The frequency-domain full-acoustic wave equation in the case of two-dimensional media with constant density can be defined as: where x denotes the horizontal orientation of the model, z denotes its vertical orientation, v(x, z) is the velocity and p(x, z, ) is the pressure field. Assuming that the medium is homogeneous, the two oneway wave equations divided from equation (1) can be expressed as: where p u and p d represent the up-going and down-going pressure wavefields, respectively, s(x, ) denotes the source wavelet, d(x, ) represents the recorded seismic data at a certain depth and Λ is the square root operator, which is regarded as the most crucial parameter for the OWDM. The typical approach to calculate wavefields in the adjacent layer using a one-way wave migration operator is p(x, z + Δz, ) = e ±iΛΔz p(x, z, ).
In this work, an eigen-decomposition method for matrices is proposed, which uses the best polynomial approximation to calculate the operator Λ. In the frequency-space domain, the Helmholtz operator L is defined as: In this equation, the diagonal elements of the matrix 2 ∕v 2 (x, z i ) consist of all velocity points of horizon coordination in the same depth, and this is why the proposed method can deal with the lateral velocity variations; the matrix 2 ∕ x 2 is a second-order finite-difference discretisation, which can be substituted into a higher-order discretisation for improving the calculation precision.
In terms of matrix transformation theory, the self-adjoint matrix L in equation (6) can be decomposed into a multiplication of eigenvalues and eigenvectors so that where the superscript T represents the matrix transpose, the matrix M denotes the eigenvalues of L, and the matrix V denotes the eigenvectors of L. When the n-order square matrix L is provided, and if there exists a set of constant { 1 , 2 , ⋯ n } ∈ M and nonzero vectors V to satisfy the equation LV = 1,2⋯n V, in that way, it can be defined M as the eigenvalue and V as the eigenvector of L.
In the framework of the matrix eigen-decomposition theory, the eigenvalues of matrix functions are functions of the eigenvalues of the original matrix and their eigenvectors remain unchanged; thus, the eigenvalues and eigenvectors of the sub-functions can be expressed as: Next, by introducing the Euler formula, the one-way wave propagator can be further expressed as: Finally, the proposed propagator for the OWDM utilising the matrix eigen-decomposition can be expressed as: (13b) Schleicher et al. (2008) compared the cross-correlation and the deconvolution imaging conditions. From their conclusions, the most precise imaging condition is the latter. In the method proposed herein, the following imaging strategy is adopted: where R(x, z, ) and S * (x, z, ) are the receiver wavefield and the conjugate of source wavefields, respectively. An issue that must be addressed is the evanescent waves in the depth wavefield extrapolation. It is well known that whether the square root operator is a positive real number depends on the operator L. When the operator L is greater than zero, significant waves propagate in the process of wavefield extrapolation. Conversely, when the square root operator becomes complex, the evanescent waves grow exponentially with the depth direction in the one-way wave migration algorithm. Evanescent waves must be suppressed during depth extrapolation because they cause instability in practice. By considering equation (10), negative eigenvalues lead to evanescent waves, and positive eigenvalues lead to physically significant propagating waves. In this study, the evanescent waves were removed by retaining the positive eigenvalues and discarding the negative ones.

Best approximation algorithm for trigonometric function calculation
Many methods can achieve a fast approximation of trigonometric functions. For instance, the search-table method works by setting the specific interval of function sampling with a fixed resolution, calculating its value, and establishing a table to be stored in memory. The corresponding address can be accessed to obtain the required trigonometric values when required. The table-checking method is efficient and straightforward, but it requires a certain amount of memory. Taylor series also can be used to evaluate the trigonometric functions. According to the Taylor theorem, a differentiable function can be approximated as the primary limit of the expansion. However, the convergence speed of the Taylor series is slow. Clearly, it is impossible to achieve high precision using a low-order approximation. Another method to approximate trigonometric functions uses least-square polynomials, which aim to minimise the square distance between the polynomial and the target function. The calculation of the LS polynomial coefficients is straightforward. Moreover, this method is easier to implement, and the approximation efficiency is much higher than that of the Taylor series. However, the LS polynomial method is not optimal as some lowerorder methods can achieve the same accuracy. In this study, the best approximation method for the cosine function is discussed and derived based on the infinite norm.
Because the function f (m) = cos(m) is continuous in the essential period [0, 2 ], the n-order polynomial approximation of the interval can be expressed as: where p n (N = 0, 1, ⋯ , N) denotes the polynomial coefficients.
To determine the polynomial coefficients, the maximum value of the difference between the polynomial f (m) and g(m) is minimised, that is (Zhou & Guo 2013): According to the numerical analysis theory, assuming that there is an extreme point m i ∈ [0, 2 ], if m = m i , it is symmetric and an additional staggered extreme point exists. Otherwise, they must be equal. In this way, the polynomial coefficients can be determined by the algorithm based on the Remez iterative method (Remez 1934).
To achieve an accuracy of 10 -3 , the order of the approximation polynomial was set to six. Considering the properties of the cosine function, the interval range of the approximating function must be limited to within [0, ∕2]. Thus, a piecewise approximation strategy was applied to ensure the most optimal approximation in the definition interval [0,2 ]. This way, the same accuracy can be achieved by a fifth order polynomial approximation, which can be written as Then, using the relation m = √ XΔz, the cosine term of the propagator in OWDM can be obtained as  (15) can be obtained by the method described previously. Finally, wavefield extrapolation can be performed by substituting equation (18) into equation (13).
As seen in figure 1, the 10th, 15th and 20th Taylor series and the nineth LS approximations are used to compare the approximation degree with a standard cosine function. It can be seen that the fifth infinite-norm polynomial approximation can achieve the same accuracy as the nineth order LS approximations and the 20th order Taylor series.
For understanding the calculation procedure of the proposed migration scheme more intuitively, the implementation process of the algorithm is as follows: Initialise the imaging data I(x, z) = 0.
(i) Loops for all depth levels.
(ii) Loops for all frequencies.
(a) Compute the operator L using equation (6). (b) Compute the eigenvalues and eigenvectors of the matrix L using equation (7). (1) For each source: Perform depth extrapolation using equation (13b) with equation (18).
GPU technology was used to speed up the algorithm evaluation and achieve high computational efficiency. Eigendecomposition is a computationally intensive part of the extrapolation process. When the wavefields at a certain depth were deduced for the initial shot gather, the adopted strategy included the storage of the eigenvalues and eigenvectors of a computed matrix using equation (7). Then, the decomposition results were applied to other shot gathers at the same depth; this way, the repeated computation of the same eigenvalues and eigenvectors of the operator L was avoided.

Impulse response testing in media with arbitrarily varying velocity
The inhomogeneous medium model measured 1000 × and classical one-way GSP scheme were used to compute the impulse responses for comparisons with the proposed eigen-decomposition scheme using the fifth infinite-norm polynomial approximation. In this experiment, the impulse response of the SFD technique was used as a typical curve, and the FWDM scheme used also served as a reference. The typical wavefronts are outlined with white dotted lines in figure 3. On a comparison, it was found both the one-way GSP impulse response curves and the proposed matrix eigen-decomposition method with the best approximation could provide almost accurate phase information when the imaging angle was less than roughly 75°. However, for greater angles, the phase information generated by the eigen-decomposition method and the used FWDM scheme matched the theoretical impulse response curve much better than the classical GSP scheme. One of the superiorities of FWDM is that there is no limitation in imaging angle. In figure 3, it can be seen that there had been no noticeable differences in imaging angles among the SFD, FWDM and proposed methods. Although the proposed eigen-decomposition method uses a low-order polynomial approximation, the highly accurate approximation of the square root operator compensates for the errors of using the minimum or average velocity as a reference velocity to approach the actual velocity. In summary, the proposed method can generate an accurate wavefield calculation of one-way propagators, especially in the case of large angles with strong velocity differences.

Impulses in media with laterally varying velocity and dipping reflection
The velocity model measured 3000 × 3000 m with a 10 × 10 m grid. A variable velocity between 1500 and 2500 m s −1 above the dipping reflector was introduced in the model. to generate a standard wavefield simulation, and the results of the classical one-way GSP algorithm and FWDM method used were compared with those obtained by the proposed scheme.
In this test, the forward simulated wavefield calculated by the SFD technique was used as a base for comparison. When the wave passes through the inclined interface, the incident and transmitted waves separate in response to the velocity difference between the two sides of the interface, as indicated by the black arrow in figure 5a. The extrapolated wavefield calculated by the migration method proposed in this study and the FWDM method used matched the reference standards. However, the extrapolated wavefields calculated by the classical one-way GSP scheme show that the incident and transmitted waves do not separate at the interface and almost stick together, as indicated by the black arrow in figure 5b. The extrapolated wavefield computed by the classical oneway GSP scheme is significantly different from that obtained from the theoretical analysis and seems to yield more imaging artifacts in the migration section. Therefore, it was concluded that the proposed migration algorithm can achieve a more precise imaging performance compared to the classical one-way algorithm.

Impulse responses in the SEG salt dome model
The SEG salt dome model has strong velocity variations in vertical and horizontal directions, which has the excellent testing ability for forward modeling and pre-stack depth migration. The description of wavefields in a medium with significant velocity changes is always the focus and main challenge of seismic inversion and imaging. An accurate description of wavefield propagation is the basis of seismic imaging and inversion; thus, it is necessary to study the wavefield calculation, and especially the description of the wavefield in complex media. In numerical experiments, the SFD technique, the classical GSP algorithm, the FWDM method used here and the migration algorithm proposed in this paper were all used to calculate the wavefield propagation in the SEG salt model. As seen in figure 6, the salt dome body velocity (4500 m s −1 ) is approximately twice that of the surrounding salt dome (2200 m s −1 ) and the velocity model shows strong velocity variations.
In this experiment, a Ricker wavelet of 30 Hz was chosen as the source and was placed at z = 0 m and x = 3000 m. The mesh interval of the velocity model was 10 × 10 m. When the seismic waves pass through the salt dome, 782 Figure 6. The salt dome velocity model. the extrapolated wavefield calculated by the classical one-way GSP scheme presents more significant errors than that calculated by the SFD technique. As indicated by the black arrows in figure 7b, the one-way GSP scheme does not accurately calculate the wavefield propagation. This is not conducive to imaging the structure below the salt dome. However, the extrapolated wavefield calculated by the algorithm in this study is consistent with the results calculated by the forward simulation of the SFD method, and closer to the FWDM method used here, as indicated by the arrow in figure 7a and 7d. Despite the more abundant waves, including the down-going waves, the turning waves and the prismatic waves generated by the FWDM method, it can be concluded that the migration method presented herein provides a more stable wavefield calculation performance in media with strong velocity variations.

Multiple shots in the Marmousi model
As is often noticed, the Marmousi model is commonly used to verify the capabilities of imaging methods. The model features steeply dipping faults with large angles, a cutting layer with sharp velocity variations and multiple anticline structures. The Marmousi model was used in this study to test the imaging performance by using multiple shots. The size of the model was 3.5 × 7.5 km, with a grid of 560 × 1200 nodes. Accurate velocity media were used to produce shot gathers by the SFD technique, as shown in figure 8a, and a smoothed velocity was used for migration, as shown in figure 8b. In this test, a wavelet with the center frequency of 60 Hz (served as a moving seismic source) was used to compute the wavefield with a sampling interval of 0.0005 s. The imaging sections generated by the classical one-way GSP scheme based on an accurate velocity model and those generated by the proposed eigen-decomposition method with the best polynomial approximation based on the infinite norm of the exact velocity model or with a smoothed velocity are plotted in figure 9a-c, respectively. No evident difference can be seen among the migration results, which means our method is not very dependent on the accuracy of the velocity model. The velocity dependence of the migration method is essential for imaging real seismic data owing to the difficulty in estimating the actual underground velocity distribution correctly. By observing the migrated results, it can be seen that all migration methods allow a good overall imaging of underground structures.
Hence, to quickly determine the differences between the two methods directly, several areas were chosen to perform a more detailed comparison, as shown in figure 9d-f. In particular, figure 9d shows the amplified part of the Marmousi velocity model as the most complex structure, and figure 9 parts e and f show the partial enlargement of the classical GSP scheme and the infinite-norm approximation migration scheme, respectively. This detailed comparison clearly shows that the proposed migration method is advantageous owing to its superior ability to deal with a wavefield at large angles and steep dipping compared with the classical oneway GSP scheme. As shown in the black dotted slanted box, the steeply dipping fault imaged by the proposed migration scheme at the contact between different strata is clearer compared to the distorted and inaccurate image generated by the classical GSP scheme. On the other hand, the region in the black dashed rectangle in the middle-deep level (a representative complex area for testing migration algorithms) remains a challenge. However, the result of the proposed migration scheme in this region appears to be more explicit and has fewer artifacts than the result of the classical GSP method. The dark arrows in the figure show that the proposed migration scheme provides more regular migration results relative to the accurate velocity model compared to the oneway GSP scheme. Overall, these results clearly demonstrate that the proposed migration scheme can indeed image complexes with large dipping angles and multi-interbedded formations.
In addition, the computational cost is a concerning problem of using the proposed scheme. As shown in Table 1, the GPU device information we used was listed, and the computation comparisons with the conventional one-way FFD propagators and the utilised FWDM scheme were also introduced. Note that for a fair comparing, the most timeconsuming part, such as the quadratic-term corrections of the FFD method, the eigen-decomposition of the proposed method and the spectral projector to suppress the evanescent waves of the FWDM method used, were all computed on the GPU graphics and their calculation results stored in the RAM.
Here, two computing strategies for implementing the proposed algorithm were compared, named strategy A and strategy B, respectively. Strategy A has been introduced already. Strategy B could be achieved by preparing all the source and receiver data for the same size in the frequency-space domain and directly extrapolating all the wavefields simultaneously. Without storing eigenvalues and eigenvectors, and less a round of loop of shot gathers, strategy B can realise a higher calculation efficiency than strategy A but with more memory requirements.

Discussions
The square root operator is a fundamental parameter in OWDM as it determines the imaging accuracy. The values of the Helmholtz operator depend on the media velocity at a given depth for a single frequency point. Therefore, existing parallel computing techniques can be implemented to improve the computational efficiency by restoring the values of the eigen-decomposition algorithm. In a three-dimensional case, the calculated quantities will sharply rise. The wavefield extrapolation operator will involve 3D tensors and their decompositions, which deserves further study.
The proposed best polynomial approximation of trigonometric functions can be applied in OWDM, as proposed in this study, but also in the depth migration scheme involved in one-way wave operators, such as the two-way depth migration (TWDM) and amplitude-preserving FWDM schemes, to improve the imaging abilities of those methods based on one-way operators. Moreover, although the migration methods based on full-wave equation have some advantages in seismic exploration, including the extensively studied RTM and the recently fast-developed TWDM, they still have some problematic or unresolved issues such as the colossal memory overheads and artifacts in RTM, the evanescent waves issue and the boundary condition problem in FWDM. That is, different migration schemes have their peculiar characteristics and applicable conditions, and no method can be 785 Where f denotes the number of frequency points, x is the horizon grids, sn denotes the total shot number, and z is the depth grids. competent for any imaging condition. Thus, as an optimised one-way wave method, with the high-angle imaging ability and lower computational costs than FWDM, the proposed scheme is still worth developing further.

Conclusions
In this study, a new scheme was introduced that integrated Euler transformation with the best polynomial approximation of the cosine function based on the infinite norm into the matrix eigen-decomposition migration framework. Subsequently, a one-way wave-equation migration scheme was deduced for imaging complex structures by using a highaccuracy one-way propagator with a low-order polynomial approximation. The reliability of the proposed scheme was demonstrated by analysing the propagation of wavefields in several models with significant lateral inhomogeneity in velocity.
In the impulse response experiments involving an inhomogeneous medium, more accurate phase information could be acquired from the proposed scheme than from the classical GSP scheme when they were compared with the theoretical curves generated by the SFD technique and FWDM method. Similarly, the wavefields transmitted from a dipping reflector in a partial variable-velocity model showed that the imaging results obtained by the proposed scheme agreed with the theoretical propagation laws. Improving the accuracy and stability of propagators in complex media is a fundamental problem. The more refined and higher-quality imaging results produced by the proposed method in the impulse response of the SEG salt model and in the multiple shot migrated results of the Marmousi model provided strong evidence of the ability and potential of the proposed scheme for solving this issue.