The method of the pseudo-depth domain is to map the velocity field from the depth domain to the vertical time domain by coordinate transformation. We use the finite difference method to solve the transformed pseudo-depth domain wave equations. The mapping relationship between the depth domain and vertical time domain is introduced in more detail below. The main method of mapping is to transform the vertical depth domain coordinates z into pseudo-depth domain coordinates τ (one-way wave vertical travel time): τ x , z = ∫ 0 z d z v s m x , z (1) Where V s m is the smooth velocity field at the corresponding vertical position. In the pseudo-depth domain, time grids are divided according to the smoothing speed of the corresponding position, which causes an uneven d τ distribution. In order to use the finite difference method to calculate, we needed to select the same size τ field sampling points artificially. In this paper, we used the spline interpolation method to obtain the pseudo-depth domain velocity field with the same vertical interval size. The interpolated velocity field was used for the following forward modeling and migration to avoid spatial aliasing; the interval Δ τ should follow the sampling condition (Alkhalifah et al., 2001). Δ τ ≤ v min 10 f max (2) Where v min is the minimum velocity in the calculation area and f max is the maximum frequency of the seismic wave. Due to coordinate transformation, the number of longitudinal sampling points is no longer equal to the number of sampling points n z in the depth domain, and the number of sampling points n τ in the pseudo-depth domain follow the condition: n τ Δ τ ≥ τ max x , z (3) Where τ max x , z is the maximum value of 1-D one-way travel time of each trace obtained by Equation 1. Only when the sampling points satisfy the conditions in Equation 3 can the data ensure integrity in the process of coordinate mapping, otherwise it will lead to missing data. Of course, after the inversion or imaging of the data, it is necessary to transform the processed data back to the depth domain. z x , τ = ∫ 0 τ v s m x , τ d τ (4)

In fact, transforming depth domain z into the pseudo-depth domain is to convert x , z into ξ , η through coordinate mapping, as shown in Figure 1. This makes it possible to use finer spatial sampling in the low-speed zone and larger spatial sampling in the high-speed zone, thereby avoiding the waste of computing resources. When each time point is calculated, the wavefield occupies the compressed grid’s number.

Next, we derived the wave equations of the VTI medium in the pseudo-depth domain. The transformation is carried out in the vertical direction as follows: ξ = x η = τ (5)

Assuming that the stress field in the Cartesian coordinate system is P z x , z , then the stress field in pseudo-depth domain becomes P τ ξ , η . The derivative transformation is written as follows: F x = ∂ F ∂ x F x y = F y x = ∂ 2 F ∂ x ∂ y (6)

By combining Equation 6 with (5), there are: ξ x = ∂ ξ ∂ x = ∂ x ∂ x = 1 ξ z = ∂ ξ ∂ z = ∂ x ∂ z = 0 η x = ∂ η ∂ x = ∂ τ ∂ x η z = ∂ η ∂ z = ∂ τ ∂ z = 1 v s m (7)

Let η x = α , Equation 7 is equivalent to the Jacobian determinant: J = ξ x ξ z η x η z = 1 ∂ τ ∂ x 0 1 v s m = 1 α 0 1 v s m (8)

From the Jacobian determinant, the diagonal elements are non-zero, which means that the coordinate τ is non-orthogonal. Therefore, the matrix tensor can be calculated: g i j = J J T = 1 ∂ τ ∂ x ∂ τ ∂ x 1 V s m 2 + ∂ τ ∂ x 2 = 1 α α 1 V s m 2 + α 2 (9) g i j = g i j − 1 = 1 det g i j = V s m 2 (10)

From the above equations, let g = v s m . The two-way wave equations can be written in the following first-order form. ∂ p ∂ t = v 2 ∇ · q ∂ q ∂ t = ∇ p (11) Where, P is the stress field and q is the strain field.

Suppose that the unit basis vector in the coordinate ξ , η system is e 1 , e 2 . When the basis vector and the matrix tensor are known, the gradient in the new coordinate system can be expressed as (Riley and Hobson, 2006; Ma and Alkhalifah, 2013): ∇ p = p ξ g 11 e 1 + p ξ g 21 e 2 + p η g 12 e 1 + p η g 22 e 2 (12)

Combining Equation 9 with Equation 12, the above equation can be simplified as: ∇ p = p ξ + α p η e 1 + α p ξ + α 2 p η + 1 v s m 2 p η e 2 (13)

For vectors q = q x e 1 + q z e 2 , divergence can be expressed as: ∇ · q = 1 g ∂ ∂ ξ g q x + 1 g ∂ ∂ η g q z (14)

Combining Equation 10 with Equation 14, the above equation can be expressed as: ∇ · q = 1 v s m ∂ ∂ ξ v s m q x + 1 v s m ∂ ∂ η v s m q z (15)

Therefore, for the acoustic isotropy medium, the first-order velocity stress equations (Virieux, 1984; Virieux, 1986) are shown as follows: 1 ρ v 2 ∂ P ∂ t = ∂ u ∂ x + ∂ w ∂ z 1 ρ ∂ P ∂ z = ∂ w ∂ t 1 ρ ∂ P ∂ x = ∂ u ∂ t (16)

P is the acoustic pressure stress field, and u and w represent the strain in the x and z directions, respectively. The form of the pseudo-depth domain derived from Equations 13, 15, 16 can be written as follows: ρ ∂ u ∂ t = ∂ p ∂ ξ + α ∂ p ∂ η ρ ∂ w ∂ t = α ∂ p ∂ ξ + α 2 + 1 V s m 2 ∂ p ∂ η ∂ p ∂ t = ρ v 2 v s m ∂ v s m u ∂ ξ + ∂ v s m w ∂ η (17)

By Taylor expansion, the discrete form is obtained: P i , j k + 1 2 = P i , j k − 1 2 + Δ t ρ v 2 v s m Δ x ∑ m = 1 L a m V i + 2 m − 1 2 , j k U i + 2 m − 1 2 , j k − V i − 2 m − 1 2 , j k U i − 2 m − 1 2 , j k + Δ t ρ v 2 v s m Δ τ ∑ m = 1 L a m V i , j + 2 m − 1 2 k W i , j + 2 m − 1 2 k − V i , j − 2 m − 1 2 k W i , j − 2 m − 1 2 k U i + 1 2 . j k = U i + 1 2 , j k − 1 + Δ t ρ Δ x ∑ m = 1 L a m P i + 1 2 + 2 m − 1 2 , j k − 1 2 − P i + 1 2 − 2 m − 1 2 , j k − 1 2 + α Δ t ρ Δ τ ∑ m = 1 L a m P i , j + 1 2 + 2 m − 1 2 k − 1 2 − P i , j + 1 2 − 2 m − 1 2 k − 1 2 W i , j + 1 2 k = W i , j + 1 2 k − 1 + α Δ t ρ Δ x ∑ m = 1 L a m P i + 1 2 + 2 m − 1 2 , j k − 1 2 − P i + 1 2 − 2 m − 1 2 , j k − 1 2 + α 2 + 1 v s m 2 Δ t ρ Δ τ ∑ m = 1 L a m P i , j + 1 2 + 2 m − 1 2 k − 1 2 − P i , j + 1 2 − 2 m − 1 2 k − 1 2 (18) Where P is acoustic pressure field, U and W are strain variables, ρ represents density, m is the finite difference order, i and j are the positions of grid points in the transverse and longitudinal directions, and a m is the finite difference coefficients.

This paper studies the wavefield characteristics and migration methods in VTI media, using the first-order velocity stress equations of VTI (Duveneck et al., 2008) under acoustic approximation: 1 ρ V P 2 ∂ P H ∂ t = 1 + 2 ε ∂ u ∂ x + 1 + 2 δ ∂ w ∂ z 1 ρ V P 2 ∂ P V ∂ t = 1 + 2 δ ∂ u ∂ x + ∂ w ∂ z ∂ u ∂ t = 1 ρ ∂ P H ∂ x ∂ w ∂ t = 1 ρ ∂ P V ∂ z (19) Where P H and P V represent the horizon and vertical stress fields, δ and ε represent the anisotropic parameters. Similar to the above method of deriving pseudo-depth domain, the first-order velocity stress equations of VTI can be deduced to pseudo-depth domain by Equations 13, 15, 19: 1 ρ V P 2 ∂ P H ∂ t = 1 + 2 ε ∂ u ∂ ξ + α ∂ u ∂ η + 1 + 2 δ 1 v s m ∂ w ∂ η 1 ρ V P 2 ∂ P V ∂ t = 1 + 2 δ ∂ u ∂ ξ + α ∂ u ∂ η + 1 v s m ∂ w ∂ η ∂ u ∂ t = 1 ρ ∂ P H ∂ ξ + α ∂ P H ∂ η ∂ w ∂ t = 1 ρ 1 v s m ∂ P V ∂ η (20)

In practice, the parameters of vertical anisotropy are quite complex, so it is necessary to combine well-logging information with petrophysical experiments to obtain different empirical equations. In the following content of model test, Gardner’s empirical equation is used to calculate the anisotropic parameter field. ρ = 0.31 ⋅ V p 0.25 ⋅ 1000 ε = 0.25 ρ 1000 − 0.3 δ = 0.125 ρ 1000 − 0.1 (21)

Similarly, for the convenience of computer expression, Taylor expansion is used to derive the discrete form: P i , j k + 1 2 = P i , j k − 1 2 + ρ v p 2 1 + 2 ε { Δ t Δ x ∑ m = 1 L a m U i + 2 m − 1 2 , j k − U i − 2 m − 1 2 , j k + α Δ t Δ τ ∑ m = 1 L a m U i + 2 m − 1 2 , j k − U i − 2 m − 1 2 , j k } + Δ t ρ v p 2 1 + 2 δ v s m Δ τ ∑ m = 1 L a m W i , j + 2 m − 1 2 k − W i , j − 2 m − 1 2 k Q i , j k + 1 2 = Q i , j k − 1 2 + ρ v p 2 1 + 2 δ { Δ t Δ x ∑ m = 1 L a m U i + 2 m − 1 2 , j k − U i − 2 m − 1 2 , j k + α Δ t Δ τ ∑ m = 1 L a m U i + 2 m − 1 2 , j k − U i − 2 m − 1 2 , j k } + Δ t ρ v p 2 v s m Δ τ ∑ m = 1 L a m W i , j + 2 m − 1 2 k − W i , j − 2 m − 1 2 k U i + 1 2 . j k = U i + 1 2 , j k − 1 + Δ t ρ Δ x ∑ m = 1 L a m P i + 1 2 + 2 m − 1 2 , j k − 1 2 − P i + 1 2 − 2 m − 1 2 , j k − 1 2 + α Δ t ρ Δ τ ∑ m = 1 L a m P i , j + 1 2 + 2 m − 1 2 k − 1 2 − P i , j + 1 2 − 2 m − 1 2 k − 1 2 W i , j + 1 2 k = W i , j + 1 2 k − 1 + α Δ t ρ v s m Δ τ ∑ m = 1 L a m Q i + 1 2 + 2 m − 1 2 , j k − 1 2 − Q i + 1 2 − 2 m − 1 2 , j k − 1 2 (22)

Claerbout (1971) proposed that the core of RTM is time consistency. The imaging conditions based on this principle of migration are given as follows: R x , y , z ; s = U x , y , z ; t ; s D x , y , z ; t ; s (23) Where, R, U, and D represent the imaging result, the wavefield of receiver, and the wavefield of source, respectively. x , y , z are the grid point position in the calculation area, s is the source, and t is the first arrival time. Since the denominator of Equation 23 is zero at the non-interface position, Claerbout (1971) used the zero delay cross-correlation of up and down traveling waves to obtain a stable imaging profile: R x , y , z ; s j = ∑ j s max ∑ i t max U x , y , z ; t i ; s j D x , y , z ; t i ; s j (24)

The RTM method is imaged by the principle of time consistency, which uses the extended source wavefield and receiver wavefield (as shown by Figure 2). The basic steps are as follows: (1) the source wavefield propagates forward; (2) the wavefield of the receiver point propagates backward; (3) zero delay cross-correlation between forward and backward wavefields are performed to obtain imaging results.

The imaging condition used in this paper is the source normalized cross-correlation imaging condition proposed by Kaelin and Guitton (2006). Using the illumination of the source wavefield to normalize the cross-correlation imaging section can effectively reduce the influence of the source on the migration results. R x , y , z ; s j = ∑ ∑ U x , y , z ; t i ; s j D x , y , z ; t i ; s j ∑ ∑ D 2 x , y , z ; t i ; s j (25)

The GPU VTI PDRTM method workflow is given in Figure 3. The specific process of VTI PDRTM model test is given here. First, the research obtained true seismic data through the VTI forward modeling of the real velocity field, using the data as the input of the inverse time extension of the wavefield at the receiver point. Second, the research used the smoothed velocity field to simulate the migration velocity field obtained by migration velocity analysis or tomography inversion. Meanwhile, the corresponding anisotropic parameters in model testing were obtained by Gardner’s empirical equation. The coordinate mapping of the obtained parameter field was used as the input of the forwarding extension of the seismic source wavefield. Third, the wave equations in the pseudo-depth domain in the VTI medium were used to carry out an inverse time extension and forward extension of the receiver wavefield and source wavefield. The extended wavefield was imaged under cross-correlation imaging conditions. Last, the final depth domain imaging profile was obtained by back-projection of the pseudo-depth domain imaging results.

GPUs have powerful parallel computing capability. CUDA-C programming technology realizes the interaction between a GPU and its central processing unit (CPU), with the serial logic judgment, data reading, and writing part being put on the CPU. Finite difference parallel computing is handed over to the GPU, which can greatly improve computing efficiency. CUDA is composed of three software levels of grid-block-thread architecture, which corresponds to the three GPU hardware levels of the device-streaming multiprocessor (SM)-core. The CUDA software levels architecture is shown in Figure 4. In this paper, we put the finite difference operation part on the GPU and transferred data between the CPU and GPU through CUDA command.

In the two-dimension case, the position of each thread is determined by the index. The index in x and z directions can be expressed as: t i d = b l o c k I d x . x * b l o c k D i m . x + t h r e a d I d x . x t i d = b l o c k I d x . y * b l o c k D i m . y + t h r e a d I d x . y

The CUDA data transmission command is as follows: c u d a M e m c p y a , b , n u m 1 * n u m 2 * s i z e o f f l o a t , c u d a M e m c p y H o s t T o D e v i c e c u d a M e m c p y b , a , n u m 1 * n u m 2 * s i z e o f f l o a t , c u d a M e m c p y D e v i c e T o H o s t

This part was based on VTI acoustic wave approximation equations to simulate in a homogeneous medium. Forward modeling was the basis of migration imaging. Therefore, only by accurately simulating the wavefield propagation process could the accuracy of imaging results be ensured. Wavefield characteristics of VTI were simulated, and the results are shown in Figure 5. It can be clearly seen that strong pseudo shear wave interference was generated in the acoustic approximate equations of VTI medium forward modeling. This interference needed to be effectively suppressed before the migration. We suppressed pseudo shear wave interference by loading the source loop near the source, assuming that the source was isotropic ε = δ = 0 or elliptical anisotropic ε = δ (Duveneck et al., 2008). This method was used to eliminate the pseudo shear wave interference in VTI imaging in the following content. The wavefield snapshot in Figure 5 adopted a model size of 6000m × 6000m, with the source set in the precise middle. The velocity and anisotropy parameters of the model were v = 4000 m s , ε = 0.316 , δ = 0.208 , respectively, and the output time was t=0.6s. It is obvious that Figure 5B has clear pseudo shear wave interference in the middle. In Figure 5C, the pseudo shear wave interference has been well suppressed. The shot record in Figure 6 adopted the method that sets shot at the surface midpoint, from which it can also be seen that the pseudo shear wave interference was suppressed.

In this part, the forward modeling records in depth domain and pseudo-depth domain were compared and analyzed. The significance of comparison in the two coordinates highlighted that the shot records systems did not change after coordinating mapping, which was an important basis for judging whether the migration method could be carried out later. Firstly, based on the real underground structure of a certain Sheng Li (SL) area, the SL velocity model (as shown in Figure 7A) was constructed, and the pseudo-depth domain transformation method was used to transform the model into the pseudo-depth domain (as shown in Figure 7B). The purpose of converting the velocity field to the pseudo-depth domain was to compare the forward modeling accuracy with the depth domain method. Secondly, in the CPU, the depth domain forward modeling and pseudo-depth domain forward modeling were carried out by using the finite difference method. The size of the velocity field was 4810m × 3060 m, with a vertical and horizontal grid spacing of 10 m. In order to facilitate comparison, the same observation system was used. Setting sources and receivers on the surface, the shot point was set at the midpoint of the ground, and there are 481 receivers on the surface with spacing of 10 m. The main frequency of the source was 20Hz, with a vertical and horizontal grid spacing of 10 m. The time sampling interval was 0.6 ms, with a maximum sampling time of 2.6s. After mapping, the vertical size changed from 306m to 201 m.

The shot records obtained by forward modeling are shown in Figure 8. Figure 8A is the results of forward modeling in the depth domain; Figure 8B is the results of forward modeling in the pseudo-depth domain. Comparing the results of the two methods, the difference is subtle. In order to compare the difference in detail, a single trace comparison of two shot records was carried out. Single trace signals at x=2000m and x=4000 m were extracted, as shown in Figures 8C, D, respectively. It can be seen that the results of the two methods coincided well with each other, except for a slight difference in amplitude. Therefore, in the later numerical test, the pseudo-depth domain method could be used for seismic wavefield imaging. In this paper, preserved shot records were used to construct the receiver point wavefield by reverse extension, and the smooth velocity field in the depth domain or the pseudo-depth domain was used as the source wavefield by forward extension. The RTM imaging conditions mentioned above were used for migration imaging.

The calculation time and memory consumption results of the single shot records obtained by the two methods in CPU were compared, which are shown in Table 1. In the process of forward modeling, the pseudo-depth domain method could improve calculation efficiency and reduce occupied memory, which saved about 22% of the calculation cost and 34% of the memory occupied separately.

This part mainly displays the imaging results of acoustic isotropy RTM, VTI RTM, and VTI PDRTM. Firstly, we established three groups of models with different characteristics for numerical simulation to verify the applicability and accuracy of imaging methods under different structures. Secondly, a set of land field data was used to reflect the application result of this method. Thirdly, according to the test time and memory usage of the model and field data, the advantages of our method were analyzed.

The smooth velocity field in Figure 9 was used to simulate the velocity field obtained by tomography or velocity analysis in practice. We used accurate forward modeling results of the velocity field as observation shot records. Among them, the calculation of vertical travel time and the transformation in pseudo-depth domain were carried out based on the smooth velocity field. The delta and epsilon parameter fields of the VTI medium were obtained by using the empirical equations mentioned above (as shown in Figures 9E–G). During the testing process, VTI parameters also needed to be transformed into the pseudo-depth domain. In the later research of this paper the VTI parameter fields were all obtained by empirical equations, and there was no difference in the structure shape with the velocity field, but only in the numerical value and physical meaning. Therefore, the related model test in VTI media mentioned in the later research of this paper no longer showed the anisotropic parameter field. For the migration of the SL model, the shot point and receiver point were set near the surface. The source adopted a 30 Hz Ricker wavelet. The target area size was 4810m×3060m, and the area size was 4810m×2010m after transforming to the pseudo-depth domain. In total, 121 shots were arranged, with a shot distance of 40 m. The receiver points adopted the mode of full receiving, with a maximum sampling time of 4s and a sampling interval of 0.6 ms. RTM imaging theory was used to SL model in the depth domain and the pseudo-depth domain, and the imaging results were displayed in the same condition. By comparing Figure 10A with Figure 10B, it can be seen that VTI RTM can depict the horizon more accurately, and the imaging results are clearer and more continuous. In the VTI medium, by comparing the imaging results of the depth domain method, the pseudo-depth domain method could achieve high-quality imaging sections and enhance the events energy of mid deep imaging. In Figures 10C, E, the x-direction and z-direction imaging profiles of the VTI medium in the pseudo-depth domain are given. After comparison, there was not much difference between them, so only the x-direction imaging profile was given in the following research of this paper..

The pseudo-depth domain method mainly performs coordinate transformation in the vertical direction. To verify the adaptability of the pseudo-depth domain method in the long offset layered structure, we built a velocity model based on marine data (as shown in Figure 11). The size of the model area was 19650m × 2010m, with a grid spacing of 10 m. The top layer was the water layer with a velocity of 2000m/s, the shot point was placed 10 m away from the surface and the receiver was set on the surface. The main frequency of the source was 30 Hz. In total, 196 shots were placed with an interval of 100 m. The maximum sampling time was 6.336s, with a sampling interval of 0.6 ms. The region transformed to the pseudo-depth domain was 19650m×1410 m (as shown in Figure 11C). Figures 11B, D are smooth velocity fields, as the input migration velocity fields were used to simulate the velocity fields processed by velocity analysis and other methods. The imaging results after migration are shown in Figure 12. From Figures 12A, C, it can be seen that the RTM results of the depth domain and pseudo-depth domain in VTI media have little difference, the imaging quality is relatively good, and the layered structure can be clearly described from the inaccurate migration velocity. By comparing Figure 12A with Figure 12B, in the imaging of long offset layered strata, imaging artifacts will exist between migration result layers if vertical anisotropy is not considered, which negatively influence the imaging result ability to converge to the true position accurately. At the position of the black arrow in Figure 11 and Figure 12, there was some structural distortion in the velocity field and imaging results after the transformation to the pseudo-depth domain. This distortion could be eliminated after transforming back to the depth domain. The numerical simulation reflected that the pseudo-depth domain method has good adaptability for seafloor data and can get good imaging results.

In addition, to verify that the pseudo-depth domain method can be applicable in more complex structural migration, and to further clarify that the inverse transformation after migration can eliminate the structural distortion caused by coordinate mapping, a more complex overthrust structure was selected as the test object. Similarly, the velocity field of the research target is shown in Figure 13. The test area was 8340m×5000m; 166 shots were placed in total, with the shot interval of 5 m and grid spacing of 10 m, and there are 834 receivers on the surface with spacing of 10 m. The maximum sampling time was 4.3s, and the sampling time interval was 0.3 ms. The results of migration are as shown in Figure 14. It can be seen from the black arrows of Figures 14A, B that the results obtained by RTM are clearer when VTI parameters are considered. The black arrow area in Figure 14B is cleaner and has a higher resolution than Figure 14A. Similar to the long offset seafloor model, in Figure 13 and Figure 14, it can be seen that structural deformation occurs in the pseudo-depth domain at the position of the yellow arrows. After the inverse transformation of the imaging results, the deformation was almost eliminated, as the yellow arrows show in Figure 14D. By comparing Figure 14B with Figure 14D, it can be found that, at the red arrows, the VTI PDRTM improves the structural imaging details at the locations where the VTI RTM is difficult to image. In complex constructions, the pseudo-depth domain method can get satisfactory imaging results.

Next, we used the land field data set to compare the results of various migration methods. There were 204 shots of land field data used for migration, 60 seismic traces for each shot, and the shot spacing and receiver spacing were both 50 m. The migration parameter field was shown in Figures 15A, C, D, with a grid number of 1056 × 510 and a space sampling interval of 12.5 m. The time sampling interval was 1 ms, with a maximum sampling duration of 6s. The source wavelet adopted a 30 Hz Ricker wavelet. The shot point and receiver point of each shot were located 12.5 m away from the surface. After the velocity field was transformed into the pseudo-depth domain, the number of grids became 1056 × 357, as shown in Figure 15B. The seismic records of 100–105 shots extracted from the shot set in which the direct wave had been cut off are shown in Figure 16. The final imaging results are shown in Figure 17. Figures 17A–C are the imaging sections of RTM, VTI RTM, and VTI PDRTM, respectively. Figure 17D shows the result of transforming back to the depth domain. It can be seen from the comparison that the migration results in Figures 17B, D are more balanced in energy, more continuous in events and clearer in structure description.

In order to make the comparison more intuitive, we selected three imaging regions at the same position for magnification in the migration results obtained by each method, as shown in Figure 18. It can be seen from Figures 18A–C, which were selected from the left red solid line box in Figure 17, that Figures 18B, C describe the horizon more clearly after considering anisotropy. Meanwhile, the pseudo-depth domain method in Figure 18C was better than the depth domain method in depicting faults (as shown by the red arrow). For the middle red solid box selected in Figure 17, the enlarged migration results also show that the migration results were more continuous in events after taking VTI parameters into account (as shown by the black arrows in Figures 18D–F). Regarding the right red solid box selected in Figure 17, at the red arrows in Figures 18G–I, the acoustic isotropic RTM could not image the undulating layer, while VTI RTM and VTI PDRTM could clearly depict the structure. In summary, after the migration processing of model and actual land data, VTI PDRTM and VTI RTM had almost the same imaging accuracy, with VTI PDRTM being better at characterizing some faults and minor structures. Then, we compared the calculation time and memory usage of various methods in respective tests.

As shown in Table 2, this paper gives the calculation efficiency and memory consumption of each region in RTM. It can be seen that VTI-RTM using GPU parallel acceleration strategy simulation could save about 70% of the computing cost compared with the CPU. Under the same conditions, the pseudo-depth domain method could reduce memory consumption by 30% and calculation time by 10–30%.