Efficient Acoustic Reverse Time Migration With an Attenuated and Reversible Random Boundary

Pre-stack reverse time migration (RTM) based on the two-way wave equation has been proved to be the most accurate seismic migration method theoretically. However, it requires reverse-order access to the wavefield calculated in forward time. In recursion computing, such out-of-order access requires that most of the recursion history should be stored on the hard disk. For massive amounts of seismic data, loading the saved wavefield data from the disk during imaging has been the bottleneck of RTM, restricting its wide application. To solve this problem, the wavefield in forward time must be reconstructed in reverse order. Although the random boundary can avoid the disk requirement by creating random velocity around the computational domain when propagate the source function. However, the random wavefield reflected from the boundary can generate unwanted artifacts in the final images. In this paper, we develop an attenuated and reversible random boundary condition which is implemented by mixing the reversible attenuation and random boundary conditions. Similar to the random boundary scheme, the proposed method just needs to save the last one or two wavefield snapshots into the memory in forward process. It then reconstructs the source wavefield in reverse order, while greatly reduces the disk input and output (I/O) requirements. Taking the attenuated property into consideration, the artificial events reflected from the boundary can be eliminated. Thus, our method can improve the imaging quality largely compared with the random boundary scheme. Numerical results demonstrate that the RTM images with our proposed attenuated and reversible random boundary condition can not only eliminate the unwanted artifacts, but also improve the computational efficiency greatly.


I. INTRODUCTION
The computational requirements of geophysical algorithms are consist of large memory size and high hard-disk access speed. With respect to the growth of seismic data and model size, the computational power of recent hardware has dramatically increased. However, the increasing in the memory size and hard disk access speed have been moderate. This relatively asymmetrical growth has prevented the widely application of some modern algorithms, such as the reverse time migration (RTM) algorithm.
Since its adoption by Baysal et al. [1], RTM quickly becomes the gold standard for high-end imaging. The core of the RTM algorithm is modeling kernel. It consists three The associate editor coordinating the review of this manuscript and approving it for publication was Wen Chen . parts: forward propagation, backward propagation, and the imaging condition. Because the source propagates from the minimum to the maximum time, while the receiver wavefield must be propagated backward in time, from the maximum to the minimum time. Ideally, the seismic wavefield should be propagated in an infinite computational domain, which is computational impossible. Instead, we must minimize the artifacts reflected from the boundaries, caused by the limited computational domain. The boundary reflections are usually attenuated by an absorbing boundary condition, such as perfectly matched layer (PML) absorbing boundary condition and damping region methods [2]- [5]. Subduing of the boundary reflections is often combined with techniques that kill the plane waves perpendicular to the computational boundary. However, all of these techniques have been proved effective for modeling seismic data but force propagation only in a single direction. For correlation imaging condition, one of the forward and backward wavefield must be stored either completely or in a check-pointed manner to the disk. The wavefield stored in the disk must then be read into memory when we correlate the forward and backward wavefield at the equivalent time positions. Although the RTM based on excitation amplitude imaging condition [6], [7] can reduce the computational time and storage requirement, the imaging time is hard to obtain.
In reality, the amount of stored wavefield data can easily reach 200 GB in the two-dimensional (2D) case and 5 TB in three-dimensional (3D) case, respectively. The wavefield of such size cannot be held in the memory of modern computers, and must be written to the hard disk and read back in when needed. The huge data input/output (I/O) cost associated with disk reading/writing (R & W) forces the Central Processing Unit (CPU) or Graphics Processing Unit (GPU) in waiting status [8], which increases the overall runtime and reduces the computational efficiency of the algorithm, especially in 3D case.
Recently, the authors of [9] listed five ways to reduce the I/O cost associated with saving huge of wavefield. For example, Symes [10] proposed checkpoint method to handle forward and backward different wavefield propagation directions, which just need to save the wavefield only at the checkpoint time. Although the checkpoint method can reduce the I/O cost significantly, the amount of the computation will be increased at least three times. Dussaud et al. [11] adopt a method which just saves the wavefield at the boundary region, which is the most popular method to reduce the I/O cost. However, in 3D RTM, the wavefield size of boundary area is still too huge to save to the disk. Yang et al. [12], [13] adopt the forward wavefield reconstruction method by interpolating significantly decimated boundary, and implement it in a staggered grid GPU. However, the saving wavefield of these methods still require significant I/O, especially for 3D RTM. Besides, references [14], [15] and [16]- [18] adopt the data compression and encoding methods to reduce the I/O cost. However, compression algorithms itself are very time-consuming.
To overcome this disadvantage, reference [19] proposed the random boundary condition to further reduce or even avoid the I/O cost associated with the wavefield. In this approach, the conventional damped region is replaced with an increasingly random velocity region. The reflected wavefield is not eliminated, but is distorted to minimize coherent correlations with the receiver wavefield. The random boundary condition achieves this in RTM imaging by making it possible to generate time-reversible wavefield that are incoherently scattered from the boundary region. The random velocity in the boundary regions is related to the velocity at the edge of the inner area. Such a boundary design effectively scatters the coherent wavefronts. In RTM with random boundary condition, the source function is first propagated from zero to the maximum time, and the wavefield of the last two time steps are stored in memory. Next, we propagate the receiver wavefield from the maximum time to zero, and reconstruct the source wavefield from maximum time to zero using the wavefield saved in the memory during the first propagation. Finally, the imaging condition can be applied on the fly and the wavefield storing is avoided. When the wavefield reaches the random boundary, the reflected wavefield is random, and has no correlation with the efficient wavefield. When we apply the crosscorrelation imaging condition, the random wavefield has little affects to the final image. Although this method can avoid the wavefield data I/O, the reflected random wavefield data has the same order of magnitude as the effective wavefield, which can introduce non-negligible artifacts in final image.
Reference [20] applied the viscoacoustic wave equation as the boundary condition in seismic modeling and inversion in time-space domain, and then implemented it in frequency domain. Zhu and Harris [21] and Zhang et al. [22] proposed their own fractional viscoacoustic wave equations with the amplitude attenuation and phase dispersion terms decoupled through different ways. Based on these equations, plenty of researches on Q-RTM algorithm and its computation efficiency have been conducted [8], [23]- [27]. Currently, the compensation accuracy and calculation efficiency can reach a satisfactory imaging result. In fact, the substance of Q-RTM is a process of wavefield reconstruction. Based on these research in Q-RTM, we adopt the lossy viscoacoustic wave equation derived by reference [23] as the boundary condition to attenuate the energy of the wavefield in the boundary area. And we also use the framework of Q-RTM to reconstruct the wavefield in forward time.
As mentioned above, we propose an attenuated and reversible random boundary condition to RTM. On the one hand, we solve the boundary area by the lossy viscoacoustic wave equation, which can attenuate the energy of the reflected wavefield from the boundary. Meanwhile, the energy of the wavefield in the inner and boundary areas can be recovered in reverse order according the framework of Q-RTM by using the wavefield at the last two time steps. On the other hand, the velocity at the boundary region is also random which is designed as the same as the random boundary condition. By mixing these two measures, the energy of the wavefield reflected from the boundary is almost one order of magnitude lower than that of the effective wavefield in the inner area. By applying the crosscorrelation imaging condition, the unwanted artifacts in RTM images using random boundary condition can be eliminated.
In this paper, all the RTM algorithms are implemented by CUDA based on GTX1080Ti with memory size of 11GB. To clarify our concept to readers, we remainder of this paper as follows: Firstly, we will review the traditional RTM with PML absorbing and random boundary conditions, and figure out its huge I/O cost. Secondly, we will introduce our new attenuated and reversible random boundary condition carefully, and implement it in acoustic RTM. Thirdly, we use a complex Marmousi model and Sigsbee2A models to show VOLUME 8, 2020 the efficiency and accuracy of our method in RTM. At last, some discussions and conclusions will be given.

II. METHOD
A. BACKGROUND RTM consists of three steps: forward propagation of the source function, backward propagation of the receiver data and application of imaging conditions. In this paper we adopt the cross-correlation imaging condition expressed as where, T represents the temporal length of the data, and x is the spatial coordinate in the x and z directions. S (x, t) is the wavefield when the source function propagates within the computation domain from t = 0 to t = T , where t is the time variable. The recorded data is injected into a second computational domain and propagated from t = T to t = 0. Then, we can obtain the backward wavefield R (x, t). By using the imaging condition given by Equation 1, the final RTM image will be obtained.

Input:
Velocity(x) of the computation domain, and Perfectly Matched Layers (PML) absorbing boundary condition for the boundary area.

Output:
The imaging result I (x). 1: for all shots do while t < T do Propagate the source function to time t = t i . and store all the wavefield S(x, t) to the disk. end while time t = T 2: Read the source wavefield S(x, t) from the disk at time t = t i . Correlate the source and receiver wavefield at constant t i using Equation 1. end while end for 3: return The imaging result I (x). Algorithm 1 show the framework of standard RTM. Note that the wavefield S (x, t) and R (x, t) are obtained in opposite temporal directions. Because the artificial boundary reflection is handled by the PML boundary condition, the wavefield of the source function cannot be reconstructed in reverse order. Therefore, S (x, t) is too large for the memory storage and must be stored to the disk. The image I (x) can be updated by Equation 1 while propagating the receiver wavefield and reading the source wavefield from the hard disk. Thus, huge I/O cost of these processes significantly reduces the computational efficiency of RTM imaging method, and become the most important barrier to limit its widely application.

B. RANDOM BOUNDARY
It has been proved that all the absorbing boundary conditions only force the propagation in a single direction, and cannot allow time reversal. Reference [19] implemented the time reversal in acoustic media by introducing a random component at the boundary of the velocity field. Because the random velocity cannot attenuate the energy of the wavefield. Thus, the wavefield can be recovered in the reverse time direction by using the wavefield of the last two time steps in source function propagation. The basic algorithm for RTM with random boundary condition can be found in Algorithm 2. Figure 1(a) shows a simple velocity model which demonstrates the construction with random boundary. The model size is 9.75 km×9.75 km with a grid interval of 4 m. The source is a Ricker function with a dominant frequency of 30 Hz, and located at the center of the model (4.875 km, 4.875 km). Figure 1 The left column in Figure 2 shows the wavefield developed in a single modeling experiment at 0.7 s (top), 1.75 s (middle) and 4.2 s (bottom), respectively. In Figure 2(c), when the Algorithm 2 Algorithm of RTM With the Random Boundary Condition Input: Velocity(x) of the computation domain.

Output:
The random boundary and imaging result I (x Reconstruct all the source wavefield t = t i using the wavefield at time t = t i . Correlate the source and receiver wavefield at constant t using Equation 1. end while end for 4: return The imaging result I (x).
wavefront of the wavefield reaches the random boundary, the wavefield become randomized and has non coherent energy. When the wavefield propagate to 4.2 s, the wavefield of the whole computational domain are all randomized. Using the cross-correlation imaging condition, the random energy reflected from the boundary will be eliminated.
In order to see if there is a coherent pattern underneath the random looking field. We conduct the experiments 50 times, each with a different random boundary. The right column of Figure 2 shows the summed results of these 50 experiments at 0.7 s (top), 1.75 s (middle) and 4.2 s (bottom). Note that averaging the results of the 50 experiments will greatly reduce the energy reflected from the boundary.
However, we can still notice the obvious noise in Figure 2(f) which will affect the final RTM image. To further analysis the remaining noises, we perform   Figure 2(f), we can notice that the wavefield is dominated by low spatial wavenumber features. But the energy still remains strong. Although the energy can be attenuated by adding the number of the sources (Figure 3(b)), it can also degrades the imaging quality of RTM. We can see the unacceptable level of low frequency noise in the images of reference [19] and the examples in this paper. Therefore, the RTM with this traditional random boundary condition cannot satisfy the requirement of high precision seismic imaging.

C. ATTENUATED AND REVERSIBLE RANDOM BOUNDARY
As mentioned above, the RTM with the random boundary condition can avoid the cost of input/output transfers to/from the disk. However, the random noise reflected from the boundary can degrade the quality of RTM. To fully exploit the advantages of random boundary that can avoid wavefield storage, while eliminating the shortcomings of strong noise that damages the quality of RTM. In this paper, we design a new attenuated and reversible random boundary condition which is shown in Figure 4. As shown in Figure 4(a), the computational domain is divided into three regions: (1) an inner area, (2) a very thin transition area that prevents the sudden reflection from (3) boundary area. The velocity of areas (2) and (3) are random which are conducted by Algorithm 2 in subsection ''random boundary''. For 2D case, areas (1) and (2) are computed by the traditional acoustic wave Equation where p = p (x, z, t) is the wavefield and v is the spatially varying velocity. In this paper, Equation 2 is solved by the finite-difference method. Area (3) where, 0 cos 2 (π γ /2), v 0 and ρ are the reference velocity and density, ω 0 is the reference angular frequency, M 0 represents the bulk modulus. As mentioned in reference [21] and [25], the second term of Equation 3 controls the amplitude attenuation of the wavefield. In area (3), the attenuation factor Q values varies from 80 to 10 with a linear variation across the boundary region shown in Figure 4 is a fractional Laplacian partial differential equation which is hard to be solved by finite difference method. Thus, we adopt the local FFT algorithm to compute these fractional Laplacian partial differential equations [28]- [30], which can implement block calculation in boundary area (3).
When the wavefield propagates from area (1) to transition area (2), its reflection is random. Additionally, when the wavefield propagates from the transition areas (2) to the boundary area (3), there will be two kinds of weak reflected wavefield. The first reflected wavefield is introduced by the abrupt change of Q, and the second one results from the different numerical calculation methods in areas (2) and (3). However, these two kinds of reflected wavefield will be randomized in transition area (2). Thus, these weak unwanted reflected wavefield will not affect the effective signal of inner area (1). And most of the energy of the wavefield will be attenuated in area (3), the energy of the wavefield reflected back to the inner area will be very weak and random. In this way, we can overcome the disadvantage of the traditional random boundary condition.
To illustrate the effectiveness of our method, we repeat the simple modeling experiment described in the section subsection ''random boundary''. The results are shown in Figure 5. Figures 5 and 2 have the same calculation parameters and display mode except the boundary condition. Comparing Figures 5(e) and 2(e), we observe that most of the reflected wavefield have been attenuated. However, the reflected wavefield with the traditional random boundary is sufficiently strong to affect the inner area. When we repeated the experiments 50 times (sources), there will be non-effective reflective wavefield which is shown in Figure 5 Figures 6  and 3, we observe that most of the reflected wavefield with our new boundary condition can be eliminated even for only a single source. Along with the increasing of imaging sources, there will be non effective reflected wavefield which is shown in Figure 6(b).
As mentioned above, during the forward propagation of RTM, the wavefield in the inner area (1) and transmit area (2) are computed by Equation 2, and the boundary area (3) is solved by lossy viscoacoustic wave Equation 3. But in backward wave propagation, when we use the cross-correlation imaging condition, we must reconstruct the wavefield in forward propagation at all time. In Q-RTM, we already can compensate the Q effect with high accuracy. Therefore, when we reconstruct the wavefield of forward propagation during backward wavefield propagation process. The areas (1) and (2)     Velocity(x) of the computation domain.

Output:
The random boundary and imaging result I (x). wave propagation using Equation 4. The whole RTM method with our new boundary condition can be found in Algorithm 3.
As is well known, when we reconstruct the wavefield of forward wavefield propagation using the wavefield at maximum time by solving Equation 4 in area (3), the process will be unstable because of the exponentially bossted high frequency ambient noise. To prevent the high-wavenumber noise from growing exponentially during the process of reconstruction, we apply a time dependent low-pass filter to the second term of Equation 4 which controls the amplitude in wavenumber domain [31]. In this paper, we select the time-variant cutoff wavenumber. The workflow of the low-pass filtering method can be summarized as follows Setp 1: Compute the reconstructed wavefield p in reverse time order (from t = T to t = 0) by solving Equation 4.

VOLUME 8, 2020
At the same time, we can obtain the corresponding wavenumber spectrum. Note, this step doesn't cause any additional computation amount because the boundary area is implemented by pseudo-spectral method.
Setp 2: Convert the wavenumber spectrum to the power spectrum. By setting a specific power threshold, we can obtain a cutoff wavenumber.
Setp 3: Take the obtained cutoff wavenumber as the regularization operator, and repeat Step 1 -Step 2, Finally, we can keep the process of wavefield reconstruction stable.
From the steps above, we can notice that our low-pass filtering method is dependent of time, the stability and accuracy will be much higher than the traditional low-pass filtering method. In addition, because only the boundary area (3) is solved by Equation 4, the increasing of computational amount can be ignored.

III. ACCURACY, STORAGE REQUIREMENT AND COMPUTATIONAL EFFICIENCY OF OUR NEW BOUNDARY A. ACCURACY ANALYSIS
In this section, we use the simple homogeneous model which has been illustrated in section ''Random Boundary'' to demonstrate the accuracy of reconstructing wavefield by using our new attenuated and reversible random boundary approach. Figure 7 shows the same traces extracted from the snapshots at 450 ms which was obtained by our new attenuated and reversible random boundary condition. The red and blue lines represent the forward and reconstructed wavefield, respectively. From which we can see that although the forward wavefield (red line) and reconstructed wavefield (blue line) cannot reach two or more orders of magnitude errors, the red and blue lines can matched well to satisfy the requirement of RTM.

B. STORAGE AMOUNT AND COMPUTATIONAL EFFICIENCY ANALYSIS
In order to shown the effectiveness of our method, we compare the storage requirement and computational efficiency with different model size between RTM with PML (RP), random (RPB) and our new boundary conditions (RRAB). The size of the models are 100 × 100, 200 × 200, 400 × 400, 800×800, 1600×1600 and 3200×3200. The boundary layers for RP, RPB and RRAB are 30, 40 and 40, respectively. There are 4000 time steps.  . The comparison of computational amount between RP, RPB and RRAB. The calculation time of PR is the summation of forward wavefield propagation, snapshots storage output (from GPU device to CPU device and then to hard disk) and input (from hard disk to CPU device and then to GPU device). And the elapsed time for RRAB (our method) is the summation of forward wavefield propagation and wavefield reconstruction. Figure 8 shows the storage requirement between RP, RPB and our new boundary conditions (RRAB). Obviously, our method can reduce the data size dramatically, which avoids frequent memory copy between GPU and CPU. In order to prove its high efficiency, we also compare the elapsed time of RP, RPB and RRAB. For the RP method shown in Figure 9, the whole elapsed time will reach 33747.76 s when the model size is 3200 × 3200, while it only takes 101.87 s for our new method.
From the comparison of speed-up ratio of RPB and RRAB shown in Figure 10, the RTM method with our new boundary, the speed-up ratio can up to 331.28. Figures 8 -10 show our method not only can avoid the storage requirement, but also improve the computational efficiency greatly.
In Figure 9, we notice that the elapsed time of our method is slightly more than that of RTM with random boundary, the reason is that the bounday area in our method is solved by pseudo-spectral method, which will spend more time than random boundary area with finite difference method. Although the speed-up ratio of our method is slightly less than that of random boundary, our method can also improve the computation efficiency greatly.

IV. MIGRATION EXAMPLES
In this section, we will use the complex Marmousi, realistic Sigebee2A models in 2D case to illustrate the accuracy of RTM using our new attenuated and revisable boundary condition.

A. MARMOUSI MODEL
The Marmousi model is used to prove that our new attenuated and reversible boundary condition can overcome the disadvantage of the traditional random boundary condition in complex media. The model size is 751 × 2301, with a discretization of 7.5 m. The velocity ranges from 1500 to 5500 m/s. The source is Ricker wavelet with a dominant frequency of 20 Hz for modeling and migration. There are 230 sources, and the recording length is 7.0 s with a temporal sampling interval of 0.7 ms. The receivers spread over the surface with a interval of 7.5 m. Figure 11(a) displays the true Marmousi velocity, Figure 11(b) stands for the smoothed migration velocity. Because the PML absorbing boundary is the best artificial boundary condition to prevent the unwanted reflected energy from the boundary, the RTM snapshots and image with PML absorbing boundary condition can be recognized as the reference one. Figure 12 displays the snapshots at 1.5 s, the source is located at the center of the surface. Figure 12(a) stands for the reference snapshot, Figures 12(b) and 12(c) show the snapshots with only one single source and average of 20 sources in forward wavefield propagation, respectively. And Figure 12(d) is the average of 20 reconstructed snapshots.  From Figure 12, we can see that the wavefield in the forward wave propagation with our new boundary condition can be reconstructed. In Figure 12(b), although the snapshot with only one single source still remains few random noises, the snapshot can reach almost the same precision as the reference one (Figure12(a)) after the average of 20 sources (Figure 12(c)). Compared with Figures 12(a) and 12(c), Figure 12(d) (average of 20 reconstructed wavefield) matches them well. Figure 13(a) presents the RTM result with the PML boundary condition (the reference image), in which the huge dataset of forward-propagated wavefield must be stored to the hard VOLUME 8, 2020 Comparing Figures 13(b) and 13(a), although the RTM image with the traditional random boundary condition ( Figure 13(b)) can obtain an imaging result without storing the huge wavefield data, strong noise reduces the signal to noise ratio (SNR) of the image. However, the image with our new attenuated and reversible random boundary condition can not only avoid the plenty of the storage, but also eliminate the unwanted noise to obtain a satisfied RTM image.
In order to make a clear comparison, we enlarged the areas enclosed by the dashed rectangles in Figure 13. The enlargements are presented in Figure 14. As indicated by the dashed arrows, the RTM image with the traditional random boundary condition (Figure 14(b)) contains strong unwanted noise that deteriorated the RTM result, especially in areas A and B. However, the RTM image obtained with our new boundary condition can accurately matched the reference image. Table 1 lists the hard disk requirements and quality evaluations of the RTM results under different boundary conditions for Marmousi model. From Table 1, although the RTM with the PML absorbing boundary condition can obtain a good result, it requires write/read operations of 69.53 GB wavefield data to/from the hard disk which reduces the computational efficiency of imaging. RTM with our new boundary can not only avoid the storage requirement, but also achieve  an acceleration of 152.3 times. Unlike with the RTM with traditional random and PML absorbing boundary conditions, the RTM with our new attenuated and reversible boundary condition can achieve a low-noise image and negligible hard disk storage of the wavefield data.

B. SIGSBEE2A MODEL
To further verify the feasibility of our algorithm to image the complex structure with high speed abnormal body, we use 2D Sigebee2A model to perform the tests. The model size is 1201 × 3201 with a spatial interval of 7.5 m. There are 320 sources evenly distributed at the depth of 7.5 m and 3201 receivers at the surface. The dominant frequency of source wavelet is 20 Hz. The record length is 9.8 s with interval of 0.7 ms and the seismic data is simulated with tenth-order spatial FD scheme. Figures 15(a) and 15(b) show the true and migration velocity models, respectively. We also conduct three RTMs, ie., RTM with PML absorbing boundary condition (Figure 16(a), referred to as the reference image), RTM with traditional In the reference image Figure 16(a), the salt and other structures can be imaged correctly with high SNR, and the interlayer and shadow layers are well clarified. However, the RTM image obtained with traditional random boundary condition (Figure 16(b)) contains the strong imaging noise, especially in the oval area, which greatly affects the resolution and SNR of the seismic imaging. In contrast, our method attenuated the wavefield in the boundary area, thereby reducing the energy of reflected wavefield (Figure 16(c)). Therefore, the RTM image obtained with our new attenuated and reversible boundary condition not only can obtain a clear image with high resolution, but also avoid the large storage problems. In order to see the advantage of our new method clearly, we also enlarged the area of Figure 16 ranged from x = 7.5 km to x = 20.25 km in horizontal direction and z = 1.5 km to z = 3.75 km in the vertical direction shown in Figure 17. In Figure 17(b), we can see the strong noise in the whole image, especially in areas A and B, which affects the quality of the image. However, the RTM image with our new method (Figure 17(c)) can obtain a good image which can match the reference image (Figure 17(a)) well. Table 2 compares the hard disk requirements and quality of RTM under different boundary condition for Sigebee2A model. Although RTM with PML absorbing boundary condition can obtain a good result, we must storage 218.75 GB wavefield to write and read from the hard disk which greatly reduces the computation efficiency and increases the calculation cost. Different from the RTM with traditional random and PML absorbing boundary conditions, the RTM with our new attenuated and reversible boundary condition can not only reduces the huge storage requirement, but also gets a speed-up ratio of 252.1 times.

V. DISCUSSION
Under the PML absorbing boundary condition, we cannot reconstruct the forward wavefield for RTM. Thus, the wavefield data must be entirely stored on the hard disk. When we use the traditional random boundary condition to the RTM algorithm, the random wavefield reflected from the boundary VOLUME 8, 2020 cannot be attenuated. Therefore, strong noise will appear in the finally images which will affect the quality of image greatly. In this paper, the random wavefield reflected from the boundary can be attenuated by a lossy viscoacoustic wave equation, when we reconstruct the forward wavefield using the wavefield at the last two time steps, the attenuated energy can be compensate by using the theory of Q-RTM in the boundary area. Thus, we can implement the advantage of RTM with traditional random boundary condition to obtain a clear image with high resolution.
However, we must clarify two key issues. Firstly, the acoustic wave equation is solved by FD algorithm in the inner area. In the boundary area, the lossy viscoacoustic wave equation is solved by pseudo spectrum method. Although the pseudo-spectrum method can increase the computation cost, the increment is negligible because the size of the boundary area is much smaller than the inner area. Thus, the additional computation amount can be ignored. Secondly, when we reconstruct the forward propagation wavefield, a low-pass filtering method should be applied to make the process of reconstruction stable. According to our test results, the process is accurate and stable in the 2D case.
The proposed method has good applicability which may be applied in other RTMs, such as Least-squares RTM and FWI in anisotropic and attenuating media, even in 3D case of these RTMs. However, the proposed method is just a beginning, and still has some problem to be addressed, for example how to maximize the randomness to minimize the reflected wavefield from boundary. We will delve further into these problems in the future.

VI. CONCLUSION
In this paper, we propose an attenuated and reversible random boundary condition. When the wavefronts of the wavefield reach the boundary area, most of the energy of the wavefield in boundary area will be attenuated due to the Q efferts. At the same time, the reflected wavefield will be very weak and random, which cannot affect the wavefield of inner area. Because the boundary area is solved by lossy viscoacoustic wave equation, we can use the theory of Q-RTM to reconstruct the wavefield of forward wave propagation.
In RTM, Firstly, we propagate the source function from time t = 0 to t = T , and the wavefield at t = T is stored in memory. Secondly, we propagate the receiver wavefield and reconstruct the forward wavefield. Then, the imaging will be obtained by cross-correlating the reciver and reconstructed wavefield. In this way, our method can avoid the huge wavefield storage requirement, and improve the computation efficiency greatly. The migration results of the synthetic data show that the RTM imaging with our attenuated and reversible random boundary condition can obtain almost the same results as RTM with the PML absorbing boundary condition in 2D case, which incurs a huge hard disk memory cost. Meanwhile, the random boundary condition is computationally efficient but introduces unwanted artifacts. Our method removes the problems and exploits the advantage of both methods, achieving high quality images and computation efficiency. In additional, our method can also be used in imaging of other areas, such as medical imaging.
i n t i p = i z * L1+ i x ; f l o a t k x _ c u t = t a p e 1 ; f l o a t k z _ c u t = t a p e 2 ; f l o a t t a p e r _ r a t i o = 0 . 2 ; f l o a t xc = k x _ c u t * L1 * dx / ( 2 * p i ) ; f l o a t z c = k z _ c u t * L2 * dz / ( 2 * p i ) ; f l o a t x s =(1− t a p e r _ r a t i o ) * xc ; f l o a t z s =(1− t a p e r _ r a t i o ) * z c ; i n t nxh=L1 / 2 ; i n t nzh =L2 / 2 ; i f ( i z >=0&&i z <=L2−1&&i x >=0&&i x <=L1−1)