Reverse Time Migration Based on the Pseudo-Space-Domain First-Order Velocity-Stress Acoustic Wave Equation

Reverse time migration (RTM) is an ideal seismic imaging method for complex structures. However, in conventional RTM based on rectangular mesh discretization, the medium interfaces are usually distorted. Besides, reflected waves generated by the two-way wave equation can cause artifacts during imaging. To overcome these problems, a high-order finite-difference (FD) scheme and stability condition for the pseudo-space-domain first-order velocity-stress acoustic wave equation were derived, and based on the staggered-grid FD scheme, the RTM of the pseudo-space-domain acoustic wave equation was implemented. Model experiments showed that the proposed RTM of the pseudo-space-domain acoustic wave equation could systematically avoid the interface distortion problem when the velocity interfaces were considered to compute the pseudo-space-domain intervals. Moreover, this method could effectively suppress the false scattering of dipping interfaces and reflections during wavefield extrapolation, thereby reducing migration artifacts on the profile and significantly improving the quality of migration imaging.


INTRODUCTION
Based on the theory of the two-way wave equation, the reverse time migration (RTM) algorithm was conceived in the early 1980s (McMechan, 1983;Whitmore, 1983). Since the wave equation does not need to be decomposed, there is no strata dip angle limitation caused by the wave equation approximation. The RTM is recognized as an ideal imaging method for complex structures and has been a popular topic in the field of geophysics (Moradpouri et al., 2017;Li et al., 2018;Zhou et al., 2018;Li et al., 2020). Chang and McMechan (1987) generalized the two-dimensional RTM to the elastic wavefield and then extended it to three dimensions McMechan, 1990, 1994). Zhang and Ning (2002) proposed multi-wave and multi-component RTM based on the eikonal equation. Sun and McMechan (2001) implemented elastic wave RTM based on the separation of Pand S-waves. Yan (2012) studied the viscoelastic tilted transversely isotropic medium wave equation RTM algorithm based on rotating staggered grids. Liu et al. (2013) achieved RTM of elastic waves in porous media based on Biot's theory. Song et al. (2015) proposed the RTM of divided-order multiples to solve the problem of imaging difficulty in the regions of low illumination based on primaries. In terms of computational efficiency and storage consumption, Liu et al. (2010) applied the graphics processing unit (GPU) for algorithmic acceleration, which greatly improved the computational efficiency of RTM. Clapp (2009) and Wang et al. (2012) used the random boundary and absorbing boundary storage strategies to reduce the consumption of storage capacity. Shi et al. (2015) analyzed the effect of random boundaries and an absorbing boundary in RTM and summarized the calculation cost and storage requirement for different boundaries and storage strategies.
After a few decades of development, the RTM technology has become increasingly mature, but it still suffers from the following problems: First, the RTM is usually achieved by using the finitedifference (FD) method with regular rectangular grids. When the underground interface model is meshed by grids, dipping interfaces and undulating surfaces only can be replaced by staircase curves, which may result in false scattering and interfacial distortion during RTM. In this respect, some scholars used variable space grids Wei, 2005, 2007;Huang and Dong, 2009), in which fine grids were adopted at regions with severe variation of medium parameters. However, this method still doesn't eliminate the limitations of rectangular grids, and it increases computational cost. Chu and Wang (2005) proposed an FD simulation method based on an irregular triangular mesh used in the finite-element method. Compared with the traditional rectangular mesh FD scheme, this method can describe undulating interfaces better, but the computational complexity increased. Besides, now there are many studies on RTM from rugged topography using curvilinear meshing or unstructured triangular meshing to get rid of the staircase approximation (Lan et al., 2014;Shragge, 2014;Liu et al., 2016;Qu et al., 2019;Liu and Zhang, 2020). Second, in the conventional RTM wavefield extrapolation based on the two-way wave equation, it produces a large number of reflection waves (back-propagating waves) at the interfaces. On the migration profile, it forms strong low-frequency noises and artifacts generated by wavepath cross-correlation with forward-and back-propagating waves, which result in low-profile imaging quality (Du et al., 2013). To reduce the influence of reflection waves, Baysal et al. (1984) deduced that the non-reflection acoustic equation can suppress the reflection waves well in the case of small incident angle under the assumption of constant impedance of the underlying medium. On the basis of the nonreflection acoustic equation proposed by Baysal, Song (2005) realized a recursive method to calculate the non-reflection scalar wave equation by introducing a wave impedance function. Willacy and Kryvohuz (2019) tried to image steep boundaries between a salt body and surrounding sediments based on the RTM using transmitted waves. He et al. (2008) developed RTM of arbitrarily wide-angle wave equations, but the imaging effect of this method is poor in shallow regions. Yoon and Marfurt (2006) introduced Poynting vector imaging conditions into RTM to realize cross-correlation of different direction wavefields, but it has a big numerical error in the regions of complex tectonics. Liu et al. (2011) proposed an imaging condition of RTM based on wavefield decomposition that separated up-going and downgoing waves by using the F-K transform; however, the method of separating wavefields required a large amount of extra calculation and storage.
The effective solution of above two problems is of great significance to improve the imaging quality of RTM. Wang et al. (2005) deduced a pseudo-space-domain scalar acoustic equation by transforming the traditional wave equation from the time-space domain to time-traveltime domain (or "traveltime domain"). This scheme not only overcomes the problem of seismic velocity interface distortion but also effectively suppresses false scattering and reflections. However, based on the second-order partial differential acoustic wave equation, Wang et al. (2005) had realized a second-order FD solution in pseudo-space domain using regular grid, which cannot meet the needs of calculation accuracy. Based on the detailed discussion of the principle of the pseudo-space-domain wave equation, this thesis derives the high-order staggered-grid FD scheme and stability condition for the pseudo-space-domain first-order velocity-stress wave equation and achieves high-precision RTM with them.

PSEUDO-SPACE-DOMAIN FIRST-ORDER VELOCITY-STRESS ACOUSTIC WAVE EQUATION
At present, the most of FD wavefield extrapolations of the acoustic wave are based on the first-order velocity-stress acoustic wave equation. It can be written as follows: where x and z represent the horizontal and vertical coordinates of the space domain, respectively, c p is the primary velocity at point (x, z), ρ is the density at point (x, z),P represents the pressure,ṽ x andṽ z represent the velocity components in the x and z directions, respectively, s(t) is the source function, with t being time. To obtain the numerical solution of Eq. 1, we usually use differences instead of differentials to approximate derivatives based on the staggered-grid technique (Vireux, 1984). The conventional FD method, which is applied to the acoustic equation, is based on rectangular grid. When the subsurface interface model is meshed, the dipping interface can only be described by using a staircase curve. It can cause false scattering in the process of wavefield extrapolation and interface distortion at the migration profile. At the same time, the two-way wave equation can generate reflected waves at the interfaces between different velocity layers. Furthermore, strong low-frequency noises and artifacts are formed on the migration profile, which lead to low profile imaging quality.
To solve above problems, a pseudo-space-domain first-order velocity-stress acoustic wave equation is proposed in this article. In RTM of acoustic wave equation, imaging about the pressureP is usually used. Therefore, under the condition that theP is not affected, assumingv x c pṽx , v z c pṽz , P P , Eq. 1 can be transformed into (2) After discretizing the continuous model into a grid model, we set the spatial unit grid length as Δξ (where ξ can represent x or z) and that the traveltime between a grid length Δξ as Δτ ξ . Then, the space grid Δξ and traveltime Δτ ξ satisfy the relationship Δξ c p Δτ ξ , where c p is the acoustic wave velocity in the grid point. The derivative of the pressure and velocity components with respect to space can be rewritten as: Then, substituting Eq. 3 into Eq. 1 yields Eq. 4, which is the pseudo-space-domain first-order velocity-stress acoustic wave equation. Usually, to solve Eq. 4 by the FD method, first we should discretize the continuous model into a grid model and then compute the "traveltime" Δτ ξ along with the grid point interval Δξ. For simplicity, Δτ ξ is called the "pseudo-space-domain interval". In the two-dimensional case, there are four pseudo-space-domain intervals at a point P(i, j), where i and jrepresent grid coordinates in the x and z directions, respectively. In the following discussion, a pseudo-space-domain interval is denoted as Δτ ± l (i, j), where l represents i or j, "−" and "+" represent the side of the smaller coordinate grid number and the side of the larger coordinate grid number, respectively. As shown in Figure 1, Δτ − i (i, j) represents the pseudo-space-domain interval between points P(i − 1, j) and P(i, j), Δτ − j (i, j) represents the pseudo-space-domain interval between points P(i, j − 1) and P(i, j), Δτ + i (i, j) represents the pseudospace-domain interval between points P(i + 1, j) and P(i, j), and Δτ + j (i, j) represents the pseudo-space-domain interval between points P(i, j + 1) and P(i, j).
Obviously, there are no velocity parameter items in Eq. 4. When the wave equation is transformed into pseudo-space domain, the original discrete space grid point velocity information is assigned to a grid line. At the same time, additional velocity information of the interfaces intersected on grid lines can be provided for computing the pseudo-space-domain intervals. Figure 2 shows a partial schematic illustration of a mesh model after the regular meshing of a velocity interface model including a dipping interface (as shown by the black solid line in Figure 2). The primary velocities at the upper and lower sides of the interface are 3,000 and 3,500 m/s, respectively, and the grid interval is 10 m. It can be seen that, after meshing the velocity interface model according to a regular rectangular grid, the dipping interface is distorted to an obvious staircase fold line. However, in the pseudospace-domain, the "propagation time" on both sides of the velocity interface is calculated according to its actual velocity and propagation distance, and the time sampling interval corresponding to the grid line is the sum of different "time of propagation" segments. Points P and Q in Figure 2 are two adjacent spatial grid points after the velocity model is divided according to a rectangular grid, and the velocity interface as shown by the black solid line intersects segment PQ at point B. In this case, the pseudo-space-domain interval between P and Q may be calculated as Δτ Δτ PB + Δτ BQ , where Δτ PB is the traveltime along with segment PB, andΔτ BQ is the traveltime along segment BQ.  Theoretically, there is no longer a distortion of the velocity interface in the pseudo-space-domain, and even the mutation of the model parameters between adjacent grid points are weakened, so it is expected that false scattering and interface reflection in the migration calculation can be reduced.
A 2Nth-Order-Accuracy Staggered-Grid FD Scheme of the Pseudo-Space-Domain First-Order Velocity-Stress Acoustic Wave Equation In the implementation of FD numerical simulations based on the pseudo-space-domain acoustic wave equation, as well as to improve the accuracy of the simulation and suppress the impact of numerical dispersion, we need to improve the accuracy of the differences. Therefore, in this study, we deduce the 2Nth-order-accuracy staggered-grid FD expression for the acoustic wave equation in the pseudo-space-domain.
In the two-dimensional case, the coordinates are denoted as (i, j) in the space-domain model discretized by sampling interval (Δx, Δz). The corresponding coordinates of the pseudo-spacedomain are (τ i , τ j ). We take zP zτ x and zvx zτ x as examples to give the propagation time interval calculation formula in pseudo-spacedomain.
where Δτ − i (i − k, j) represents the pseudo-space-domain interval between points (i − 1 − k, j) and (i − k, j), and Δτ −(m−1/2) i represents the propagation time interval between (τ i , τ j ) and the grid point where v x is located. Similarly, the pseudo-space-domain propagation time interval of P and v z in the direction τ z can be calculated separately in a similar manner as described above.
Using the propagation time interval shown in Eqs 5, 6, the 2Nth-order-accuracy expansion of the first-order derivative of the P with respect to the variable τ x can be obtained. The P in the pseudo-space-domain is assumed to have a 2Nth-order derivative.
where P (n) represents the nth-order derivative of P. The above 2N equations are multiplied by c x 1 , c x 2 , . . ., c x 2N−1 , c x 2N , respectively, and then added and simplified as To resolve the first-order derivative FD scheme of P at τ x τ i+1/2 , Eq. 11 needs to satisfy the algebraic relationship that the coefficient of the first derivative is one and the other derivative is 0 except at the first order. Therefore, according to the coefficient relationship between derivatives, we can get the FD coefficients c x m (m 1, 2, . . . , 2N−1, 2N). and c x 2k−1 ≠ − c x 2k (k 1, 2, . . . , N−1, N). From the above analysis, it can be seen that the FD coefficients of the numerical simulation in the pseudo-space-domain are related to the grid velocity and the size of the grid. Even with the same difference order, the FD coefficients are different corresponding to various grid velocity and sizes.
By substituting the FD coefficients into Eq. 11, a 2Nth-order difference expression for the first derivative of P at τ x τ i+1/2 can be written as: Similarly, we have a 2Nth-order difference expression and FD coefficients c z m (m 1, 2, . . . , 2N−1, 2N) for the first-order derivative of P at τ z τ j+1/2 , a 2Nth-order difference expression, and FD coefficients c x vm (m 1, 2, . . . , 2N−1, 2N) for the first-order derivative of v x at τ x τ i , and a 2Nth-order difference expression and FD coefficients c z vm (m 1, 2, . . . , 2N−1, 2N) for the first-order derivative of v z at τ z τ j .
Substituting the above difference expressions for P, v x , and v z at τ x or τ z and the second-order difference expressions for the first-order derivative of P, v x , and v z at time into Eq. 4, we can obtain where k represents discrete time points, which satisfy t kΔt (where Δt represents a discrete time interval). In the following, a homogeneous model is used to show the effect of the pseudo-space-domain high-order FD scheme on dispersion suppression. The size of the model is 2000 × 2000 m. The spatial sampling interval is 10 × 10 m, and the primary wave velocity is 2,500 m/s. The source location is 1,000 m, 1,000 m. As a source wavelet, we adopt the Ricker wavelet whose dominant frequency is 35 Hz, and the time sampling interval is 0.25 ms. Snapshots based on simulations of the pseudo-space-domain with different orders of FD operator are shown in Figure 3.
It can be seen from the snapshot shown in Figure 3 that, for the pseudo-space-domain FD scheme with different orders of FD operator, when the spatial grid interval, model velocity, and wavelet dominant frequency are the same, the higher the order of FD operator is, the weaker the dispersion is.

Stability Condition for the Pseudo-Space-Domain First-Order Velocity-Stress Acoustic Wave Equation
First, we define the pseudo-space-domain plane harmonic variables u u u 0 e iωnΔt e ikτ x jΔτx e ikτ z kΔτz , where u 0 represents the initial wavefield, ω represents the circular frequency, k τx and k τz represent wave numbers in the τ x and τ z directions, respectively; n, j, and k represent coordinates of discrete grid points in the t, τ x , and τ z directions, respectively; Δτ x and Δτ z represent pseudo-space-domain intervals in the τ x and τ z directions, respectively; e stands for the base of the natural logarithms, and i represents the imaginary unit in this section.
According to the above equation, one can get the following relationships: Substituting the above formulas into the first-derivative difference expression gives Furthermore, the second-derivative expression for τ x can be obtained as According to the definition of the propagation time in the pseudo-space-domain and the sampling theorem, when the maximum wave number for τ x is obtained, the following relationships hold: Δτ +m x k τx (m − 1/2)π and Δτ −m x k τx (m − 1/2)π. Therefore, the above equation can be converted into Similarly, the second-derivative expression for τ z can be obtained as and the second-derivative expression for t can be obtained as Because the left side of the equation above satisfies 0 ≤ sin 2 ω Δt 2 ≤ 1, and under the assumption that the differential coefficients in the τ x and τ z directions are equal, the following relation holds: where c m c x m c z m (m 1, 2, . . . , 2N − 1, 2N). Frontiers in Earth Science | www.frontiersin.org October 2021 | Volume 9 | Article 690513 6

Perfectly Matched Layer Boundary Conditions of the Pseudo-Space-Domain First-Order Velocity-Stress Acoustic Wave Equation
In the central wavefield calculation region, FD numerical simulation of the pseudo-space-domain first-order velocitystress acoustic wave equation with 2Nth order in pseudo-space and second-order in time can be realized by applying Eq. 13. In the artificial boundary region, to effectively suppress the artificial boundary reflection, absorbing boundary processing is required. The perfectly matched layer (PML) boundary conditions for the first-order velocity-stress acoustic wave equation in the pseudospace-domain are given below.
Because the PML attenuation term is independent of the partial derivative of wave equation, the space domain partial derivative in the equation is transformed into a pseudo-spacedomain partial derivative; meanwhile, the space domain attenuation factors d x and d z are transformed into pseudospace-domain attenuation factors d τx and d τz . The PML boundary conditions for the pseudo-space-domain acoustic wave equation are as follows: where P x and P z represent the components of P in the τ x and τ z directions. d τ x and d τ z represent the attenuation factors in the τ x and τ z directions, which are given by where τ m represents the normal pseudo-space-domain propagation time interval from the point in the PML layer to the outer edge of the center wavefield, R represents the theoretical reflection coefficient for the PML layer (ranging from 10 -5 to 10 -7 ), and τ L represents the pseudo-space-domain PML layer thickness. As is shown in Figure 4, when wavefield calculations are performed, d τ x 0, d τ z 0 in the center wavefield area; d τ x 0, d τ z d(τ m ) in PML areas one and PML area 4; d τx d(τ m ), d τz 0 in PML areas two and PML area three; and d τx d(τ m ), d τz d(τ m ) in the corner area. We can write the attenuation factors d τx v x , d τz v z , d τx P x , and d τz P z in the PML boundary conditions into a differential form, and by substituting them with the difference expression of each first-order derivative into Eq. 23, we can derive  A uniform medium model is used to verify the effectiveness of the PML boundary conditions of firstorder velocity-stress acoustic wave equations in the pseudo-space-domain for eliminating artificial boundary reflections. The horizontal and vertical lengths of the model are 2000 and 2000 m, respectively, and the grid interval is 5 m. The primary wave velocity in the model is 2,500 m/s, and the density is 2000 kg/m 3 . The source location is (1,000 m, 1,000 m), and the source wavelet employs Ricker wavelet with a dominant frequency of 35 Hz. To avoid the effects of numerical dispersion, the pseudo-space-domain FD order is 10th order. The snapshot of the wavefield extrapolation process at 0.82 s is shown in Figure 5, where Figure 5A shows the snapshot of the left boundary without the PML boundary conditions, and Figure 5B shows a snapshot of the left boundary with PML boundary conditions.
To further illustrate the boundary absorption effect of the PML boundary conditions in the pseudo-space-domain acoustic wave equation, the left boundary reflection wave corresponding to a depth of 1 km in the wavefield shown in Figure 5 is magnified and displayed and is compared with the conventional acoustic wave equation wave based on the same simulation parameters. As shown in Figure 6, the black solid line is the left boundary reflection wave absorbed by the PML boundary condition of the pseudo-space-domain acoustic wave equation, the red dotted line is the left boundary reflection wave absorbed by the PML boundary condition of conventional method, and the blue dotted line is the reflected wave of left boundary without attenuation by PML boundary conditions. It can be seen that the amplitude of the boundary reflection wave after absorption by the pseudospace-domain PML boundary condition is basically the same as that obtained by the conventional PML boundary condition, and it is obviously weaker than the amplitude of the uncompressed boundary reflection wave.

Normalized Cross-Correlation Imaging Conditions
In this study, normalized cross-correlation imaging conditions (Kaelin and Guitton, 2006) are used in the RTM. The realization process employs the zero-delay cross-correlation of the source wavefield to normalize the zero-delay cross-correlation between the forward time wavefield and the reverse time wavefield as where T is the total recording time.
(v) F is the forward time wavefield, and (v) R is the reverse time wavefield. When using the above imaging conditions, it is usually necessary to save the forward time wavefield at each time. However, when all the wavefield values are stored on the storage medium, large amount of memory storage space is required and also a long data access time. To overcome this problem, in this study, we implement an effective boundary storage strategy (Clapp, 2008;Wang et al., 2012) based on PML boundary conditions in the pseudo-space-domain. This entails storing the wavefield value of the N-layer grid point (the FD order is 2Nth order) that is adjacent to the central wavefield on each PML boundary during the forward time source wavefield extrapolation, as well as the central wavefield value at the last moment. These boundary wavefield values are taken out as a boundary condition when extrapolating along reverse time, and then, the source wavefield can be rebuilt in the time iteration. Although this method requires the forward time source wavefield extrapolation in advance, this can effectively reduce the storage requirements for the RTM in pseudo-space-domain. The main purpose of this experiment is to test the validity of pseudo-space-domain acoustic wave equation RTM in solving velocity interface distortion and suppressing false scattering and reflected waves.
The experiment used a two-layer velocity model with a dipping interface, as shown in Figure 7A. The horizontal and vertical lengths were 4,000 m and 2000 m, respectively. The primary wave velocities at the upper and lower sides of the interface were 2,500 m/s and 3,500 m/s, respectively. The density was 2000 kg/m 3 . The grid model obtained by meshing this interface model with vertical and horizontal grid intervals of 10 m is shown in Figure 7B. It can be seen that the original smooth velocity interface has become an obviously stepped interface (white arrow in the figure). In the experiment, a geometry with a fixed position of receivers and changeable source position was established. The shot point was between 500 and 3,480 m, the interval between the shots was 20 m, and a total of 150 shots was made. There were 401 receivers per shot, and each receiver was located between 0 and 4,000 m. The interval between receivers was 10 m. The depths of shots and receivers were both 10 m.
Obviously, to verify the effectiveness of a migration method in solving velocity interface distortion, it is necessary to ensure that the acquired seismic record is accurate. The shot records required for this experiment were obtained by using FD wave equation forward modeling. In theory, only by using a small enough grid spacing can ensure that the obtained shot records are relatively accurate. Therefore, in this study, we first used a model with 1 m grid intervals in both vertical and horizontal directions for forward modeling. (Note that, even if the number of grid points is only doubled, this can cause a huge increase in the amount of   computation. Therefore, regardless of whether one realizes forward modeling in actual processing or migration and inversion, it is generally unrealistic to use such a small grid interval.) The simulation used a Ricker wavelet with a dominant frequency of 35 Hz. The FD order was 16th order in space and second-order in time, and a total of 150 shots of seismic records was obtained. The record of the 76th shot is shown in Figure 8A. For comparison, a record obtained with the grid interval of 10 m is shown in Figure 8B, which existed strong artificial scattered waves. Based on the model of 10 m grid interval as shown in Figure 7B, the FD algorithm for the conventional and pseudospace-domain acoustic wave equations is used for RTM with second-order accuracy in time and sixteenth-order accuracy in both space and pseudo-space. (Note that the pseudo-spacedomain RTM needs to add velocity interface information to calculate the traveltime between two grid points.) The migration profiles are, respectively, shown in Figures 9A,B.
To more intuitively compare the morphology of the dipping interface in the migration profile of the two methods, the event in the elliptical region in Figures 9A,B is magnified, as shown in Figures 10A,B. It can be seen that the shape of the dipping interface in the RTM profile of the conventional acoustic wave equation (the red dotted line in Figure 10A) is significantly distorted compared with the real interface morphology (the red solid line in Figure 10A). The interface shape in the RTM profile of the pseudo-space-domain acoustic wave equation is basically consistent with the real interface morphology (the red solid line in Figure 10B). This demonstrates that pseudo-spacedomain acoustic wave equation RTM can effectively solve the distortion problem of the velocity interface. Figures 11A,B show a snapshot of the forward time wavefield of the 76th shot at 0.9 s. It can be seen that there are obvious false scatterings in the wavefield of the conventional acoustic wave equation (as shown in the elliptical region in Figure 11A), and there are no obvious false scatterings in the wavefield of the pseudo-space-domain acoustic wave equation (as shown in the elliptical region in Figure 11B). By comparing the interface reflections at the arrows in Figure 11, we can see that the

Reverse Time Migration of the Marmousi Model
The Marmousi model (shown in Figure 12) is a grid velocity model of complex tectonic with numerous velocity interfaces, steep dip structures, and dramatic velocity changes. The model size is 9,200 m * 3,000 m, respectively. The horizontal and vertical grid spacings are, respectively, 5 and 4 m. In the experiment, the unilateral shot geometry used was a seismic source located at the right side and receivers located at the left side. The interval between the shots was 25 m, with 426 shots in total. There were 104 receivers per shot. The depths of the shot and the receivers were both 8 m. The source wavelet used a Ricker wavelet with a dominant frequency of 35 Hz. Synthetic seismograms were simulated by the conventional acoustic wave equation FD method whose FD order was second-order in time and eighth-order in space.  Frontiers in Earth Science | www.frontiersin.org October 2021 | Volume 9 | Article 690513 11 We applied the conventional acoustic wave equation and pseudospace-domain acoustic wave equation FD (with FD order being second-order in time and eighth-order in space and pseudo-space) to perform RTM. Figures 13A,B show the snapshots of the 138th shot based on the two methods at 1.9 s. Figures 14A,B show the RTM profiles based on the two methods, respectively. Figure 13 shows that the pseudo-space-domain acoustic wave equation has a significantly weakened reflection wave in the wavefield compared with the conventional acoustic wave equation. It indicates the validity of the pseudo-space-domain acoustic wave equation in suppressing reflected waves during wavefield extrapolation.
Comparing the local migration profiles in the red rectangular region in Figure 14, we can see that the pseudo-space-domain acoustic wave equation has a clearer structure and better continuity of the event where the arrows point the conventional acoustic wave equation. This demonstrates that the imaging quality of RTM by using the pseudo-space-domain acoustic wave equation is better than that obtained by using a conventional acoustic wave equation.

CONCLUSION AND DISCUSSION
Based on the first-order velocity-stress acoustic wave equation in the pseudo-space-domain, we derived a 2Nth-order-accuracy staggered-grid FD expression and its PML boundary condition, deduced the stability conditions of the staggered-grid FD expression, and realized RTM in the pseudo-space-domain. At the same time, numerical experiments were carried out based on a dipping interface model and the Marmousi model. Experimental results were as follows: 1) The RTM method based on the pseudo-space-domain acoustic wave equation could solve the problem of velocity interface distortion that appears in the conventional RTM profile. 2) Wavefield extrapolation based on the pseudo-spacedomain acoustic wave equation could significantly weaken interface false scattering and reflection waves, thereby further improving the quality of the migration imaging.
Of course, the high-order FD RTM method for the acoustic wave equations in the pseudo-space-domain is not ideal for reflection waves suppression under non-normal incidence. Therefore, the focus of future research work will be the further improvement of the reflection waves suppression effect of the method and accuracy of RTM imaging, along with developing it into the elastic wave equation and the RTM of the threedimensional wave equation. Frontiers in Earth Science | www.frontiersin.org October 2021 | Volume 9 | Article 690513