Depth-Extrapolation-Based True-Amplitude Full-Wave-Equation Migration from Topography

: The lack of an initial condition is one of the major challenges in full-wave-equation depth extrapolation. This initial condition is the vertical partial derivative of the surface waveﬁeld and cannot be provided by the conventional seismic acquisition system. The traditional solution is to use the waveﬁeld value of the surface to calculate the vertical partial derivative by assuming that the surface velocity is constant. However, for seismic exploration on land, the surface velocity is often not uniform. To solve this problem, we propose a new method for calculating the vertical partial derivative from the surface waveﬁeld without making any assumptions about the surface conditions. Based on the calculated derivative, we implemented a depth-extrapolation-based full-wave-equation migration from topography using the direct downward continuation. We tested the imaging performance of our proposed method with several experiments. The results of the Marmousi model experiment show that our proposed method is superior to the conventional reverse time migration (RTM) algorithm in terms of imaging accuracy and amplitude-preserving performance at medium and deep depths. In the Canadian Foothills model experiment, we proved that our method can still accurately image complex structures and maintain amplitude under topographic scenario.


Introduction
Seismic exploration is a fast and intuitive exploration tool before drilling. Almost all oil companies rely on seismic exploration results to determine the location of exploration and development wells. In recent decades, with the development of digital technology, seismic exploration can provide more and more information.
Seismic migration is an important part of seismic data processing, which can transform seismic signals into images to intuitively show underground structural information. Moreover, imaging results of true-amplitude migration can provide lithological information. The true amplitude values of the reflectors are of great benefit to seismic attribute inversion, such as AVO (amplitude-versus-offset) inversion.
Seismic exploration on land, faces challenges from complex surface conditions, such as severe lateral velocity variation and complex terrain. Seismic imaging in this complex situation has higher requirements for migration algorithms. Among the existing migration algorithms, the ray-based methods [1][2][3] are the most widely used. These methods have good adaptability to the topographic surface, and the computation cost is low, but ray-based methods have limitations for complex structure imaging. The one-way wave equation migration [4][5][6] performs better than the ray-based method in complex structure imaging. However, the one-way wave equation method is less adaptable to the topographic surface, and it is not as good as the two-way wave equation migration method in terms of imaging accuracy. As a method based on the full wave equation, the reverse time migration (RTM) method is more advantageous in seismic imaging under topographically influenced scenarios [7]. The full-wave-equation migration method not only has advantages in structural imaging, but also has better performance in terms of amplitude preservation. Because of these advantages, the RTM method has developed rapidly and has become the research focus of scholars. Many RTM from topography strategies [7][8][9] have been proposed.
As we all know, the basis of the migration method is to solve the wave equation. RTM solves the full wave equation along the time axis. Because the wavefield extrapolates along the time axis, the calculated wavefield contains multi-scattering waves. This is beneficial for RTM imaging complex structures, but it may also bring more crosstalk artifacts, such as false events under a salt dome caused by internal multiple reflections [10]. Moreover, the RTM method assumes that the wavefield value outside the maximum sampling time is zero [11]. This assumption can affect the amplitude-preserving performance of the method to a certain extent.
Unlike RTM algorithms, depth-extrapolation-based full-wave-equation migration methods solve the full wave equation in depth. In 1983, Kosloff and Baysal [12] first introduced a full-wave-equation downward continuation scheme. They used a low-pass filter in the wavenumber domain to suppress evanescent wave energy. However, this strategy loses part of the useful waves, resulting in poor imaging for steep reflectors. Sandberg and Beylkin [13] proposed a full-wave-equation depth extrapolation method, which can suppress evanescent waves by using a more advanced method of spectral projection. You et al. [14] analyzed the weaknesses of the existing methods and introduced a full-wave-equation depth extrapolation scheme based on the dual-sensor acquisition system. In their experiments, they demonstrated that the full-wave-equation migration method based on depth extrapolation has better amplitude-preserving performance than that of the conventional Kirchhoff and RTM migration methods. You and Cao [10] derived a depth-extrapolation-based full-wave-equation migration using matrix multiplication, and applied parallel computing to the core operator to improve computational efficiency.
Theoretically the full-wave-equation migration based on depth extrapolation is similar to the imaging accuracy of RTM because they are all derived from the full-wave equation. Moreover, compared with RTM, depth-extrapolation-based full wave equation migration has some potential advantages, such as less storage requirements and possible cleaner imaging results [10].
However, there are some difficulties in implementing the full-wave-equation depth extrapolation. Solving the full wave equation in depth requires two initial conditions: the wavefield value of the surface and its vertical partial derivative. The traditional seismic acquisition system can only provide the wavefield value of the surface. The lack of another initial condition is the first problem faced by full-wave-equation depth extrapolation. To solve this problem, scholars [12,13,15] estimated the vertical partial derivative of the wavefield from the wavefield value by assuming constant-velocity near the surface. On land, this assumption is often not satisfied. To move beyond this assumption, You et al. [14] proposed a dual-sensor seismic acquisition system. This acquisition system collects two sets of wavefield data at different depths, so the vertical partial derivative of the wavefield at the surface can be provided. However, it is very difficult to use such a complex acquisition system in complex terrain. Moreover, the economic cost of this acquisition system is very high. These solutions that provide the vertical partial derivative, limit the application scenarios of the migration methods based on full-wave-equation depth extrapolation.
In this paper, we propose a new vertical partial derivative calculation method. This new method is not limited by the topographic surface with lateral velocity variation, and makes it possible to apply the full-wave-equation depth extrapolation in irregular topography scenarios. Based on this, we implement a depth-extrapolation-based full-waveequation migration from topography by using direct downward continuation. Numerical experiments prove the structural imaging performance and the amplitude-preserving performance of our proposed method.

Methods
There are two challenges in full-wave-equation depth extrapolation. First, solving the full wave equation in the depth direction requires the vertical partial derivative as an initial condition, but the conventional seismic acquisition system only receives the surface wavefield value and cannot provide the vertical partial derivative. Second, we need to eliminate the evanescent wave generated in the continuation process.

Estimation of the Vertical Partial Derivative of the Surface Wavefield
Among the existing methods, the most used method to provide the vertical partial derivative has been proposed by Kosloff and Baysal [12]. They assume that the velocity near the surface is laterally uniform. Under this assumption, we can calculate the vertical partial derivative from the wavefield value of the surface. In reality, the surface conditions on land are usually more complex, such as severe lateral velocity variation and complex terrain. In order to achieve full-wave-equation depth extrapolation under complex topography, we propose a new method to calculate the vertical partial derivative of the surface wavefield.
The energy of the seismic wave is emitted from the local source, and the wavefront always has a certain degree of curvature. However, at a sufficient distance from the source, the wavefront is flat enough so that the plane wave approximation is locally correct.
The seismic wave starts from the source and reaches the ground after being reflected by the reflectors. Figure 1 shows a microscopic part of the surface: point A is on the surface, point B is very close to point A. The distance between two points is r. The angle between the raypath and the horizontal surface is θ. From the plane wave equation we can get the following: where U is the wavefield, k = w/V, and V is the P-wave velocity. From Equation (1)  Numerical experiments prove the structural imaging performance and the amplitude-preserving performance of our proposed method.

Methods
There are two challenges in full-wave-equation depth extrapolation. First, solving the full wave equation in the depth direction requires the vertical partial derivative as an initial condition, but the conventional seismic acquisition system only receives the surface wavefield value and cannot provide the vertical partial derivative. Second, we need to eliminate the evanescent wave generated in the continuation process.

Estimation of the Vertical Partial Derivative of the Surface Wavefield
Among the existing methods, the most used method to provide the vertical partial derivative has been proposed by Kosloff and Baysal [12]. They assume that the velocity near the surface is laterally uniform. Under this assumption, we can calculate the vertical partial derivative from the wavefield value of the surface. In reality, the surface conditions on land are usually more complex, such as severe lateral velocity variation and complex terrain. In order to achieve full-wave-equation depth extrapolation under complex topography, we propose a new method to calculate the vertical partial derivative of the surface wavefield.
The energy of the seismic wave is emitted from the local source, and the wavefront always has a certain degree of curvature. However, at a sufficient distance from the source, the wavefront is flat enough so that the plane wave approximation is locally correct.
The seismic wave starts from the source and reaches the ground after being reflected by the reflectors. Figure 1 shows a microscopic part of the surface: point A is on the surface, point B is very close to point A. The distance between two points is r . The angle between the raypath and the horizontal surface isθ . From the plane wave equation we can get the following:  The wavefield gradient can be obtained from Equation (2), The wavefield gradient can be obtained from Equation (2), When r approaches 0, we have Since where X and Z are the horizontal and vertical coordinates. We have ∂r ∂z Substituting Equation (6) into Equation (4), we obtain In the field, geophones for reflected wave seismic exploration are generally arranged near the source. So the distance between the source and the geophones can be considered small compared to the distance from the source to the elastic interface which usually reaches up to several kilometers. The angle between the reflected wave and the horizontal surface is close to 90 • . Under this assumption, θ = π/2, Equation (7) can be written as According to Equation (8), the vertical partial derivative is calculated locally and trace by trace, hence it can handle surfaces with topography and lateral velocity variations.

True-Amplitude Full-Wave-Equation Depth Extrapolation for Migration
For full-wave-equation depth extrapolation, eliminating the energy of the evanescent wave is the focus of the downward continuation process, because the exponentially increasing evanescent wave energy can make the numerical solution blow up. After researching the existing solutions, we use the projector PLP [13] to remove it. In the space-frequency domain, we have where d(x, z = 0, t) is the recorded seismic data at z = 0 m, FT means Fourier transform, u z (x, z = 0, ω) is the vertical partial derivative calculated by Equation (8), P = (I − sign(L))/2 (for details of operator P, see [16]).

Direct Downward Continuation from Topography
Direct downward continuation is a favorable method to solve the problem of irregular surface. The process is simple and it is also applicable to the situation of large fluctuations in the surface elevation. Ref. [4] first proposed the concept of wavefield downward continuation based on accumulating step by step. According to this concept, we achieve the full-wave-equation depth extrapolation for migration by using direct downward continuation from topography. The specific steps are as follows: (1) First, mesh the model, select the height of the highest point on the surface to establish a datum; then use the surface velocity to fill the space between the datum and the surface; and set the initial value of the source and receiver wavefields to zero.
(2) Starting from the datum, at each depth, first check whether the receiver or source point exists, if so, add the recorded wavefield or source wavelet to the corresponding wavefield; then use Equation (9) to downward continue the source and receiver wavefield.
(3) Filter the wavefield calculated at each depth to eliminate the wavefield energy above the surface [4].
(4) Update the image using the cross-correlation imaging condition. Our method uses common shot data. We perform the above steps for each shot, and superimpose the results of all shots to get the final imaging result.

Examples
In this part, we use several experiments to test the imaging performance of our proposed method.

Horizontal Layer Model
In order to achieve the full-wave-equation depth extrapolation for migration, we have three methods to provide the vertical partial derivative of the surface wavefield. The first is the traditional method of using the surface wavefield to calculate the vertical partial derivative under the premise that the surface velocity is uniform; the second is to use the double-layer wavefield recorded by the dual-sensor acquisition system to calculate the vertical partial derivative; the third is the method we proposed.
In this experiment, three migration methods are used to compare the imaging results, and we call them Method 1, Method 2, and Method 3. The only difference between these three methods is that they use different ways to provide the vertical partial derivative of the surface wavefield. Method 1 is to use the traditional method to calculate the derivative, Method 2 is to use the dual-sensor acquisition system to provide the derivative, and Method 3 is our proposed depth-extrapolation-based full-wave-equation migration method which uses Equation (8)  For each shot, the first receivers are put at x r j,i = x s i + jdx, j = 0, 1, . . . , 119 on the surface, the second receivers are placed 5 m directly below the first receivers (Method 2 uses a dual-sensor acquisition system, so we set up the second receivers). Figure 2b-d shows the imaging results obtained using three different methods. In order to study the amplitudepreserving performance, we extract the amplitude from the target reflector in the imaging result to estimate the reflection coefficient and compare the calculated reflection coefficient with its corresponding theoretical reflection coefficient. We plot the calculated reflection coefficient curves of the three imaging results and their corresponding theoretical coefficient curves in Figure 3a,c,e. Figure 3b,d,f show the relative error between the calculated reflection coefficient and the theoretical reflection coefficient.
Appl. Sci. 2021, 11, 3010 6 of 14 second receivers are placed 5 m directly below the first receivers (Method 2 uses a dual-sensor acquisition system, so we set up the second receivers). Figure 2b-d shows the imaging results obtained using three different methods. In order to study the amplitude-preserving performance, we extract the amplitude from the target reflector in the imaging result to estimate the reflection coefficient and compare the calculated reflection coefficient with its corresponding theoretical reflection coefficient. We plot the calculated reflection coefficient curves of the three imaging results and their corresponding theoretical coefficient curves in Figure 3a,c,e. Figure 3b,d,f show the relative error between the calculated reflection coefficient and the theoretical reflection coefficient. From Figure 2, we can see that the imaging effects of the three methods are very close, and the position of the reflective layer is accurate. In Figure 3, we find that the calculated reflection coefficient curves of Method 2 and 3 are in good agreement with their corresponding theoretical curves, but the calculated reflection coefficient curve of Method 1 is inconsistent with the theoretical curve. Combining Figure 3 and Table 1, we can quantitatively analyze that the relative error of Method 1 is larger than that of the other two methods, and the performance of Method 2 and Method 3 are relatively close. These results indicate that Method 2 and 3 have better amplitude-preserving performance in this experiment.
In this experiment, there is a lateral velocity change in the upper layer of the velocity model. The vertical partial derivative calculation method used in Method 1 needs to assume a uniform surface velocity. When this assumption is not satisfied, there will be a non-negligible error in the calculated derivative. It will affect the amplitude-preserving performance of the imaging result. The performance of Method 2 and Method 3 are very close, and both methods can handle lateral surface velocity variation. However, these two methods are based on different acquisition systems: Method 2 is based on a dual-sensor acquisition system, and Method 3 is based on a conventional acquisition system.
All in all, compared with traditional methods, our proposed vertical partial derivative calculation method can handle lateral surface velocity variation. Compared with the method based on a dual-sensor acquisition system, our imaging effect is similar, but our acquisition system is simpler and more practical.

Marmousi Model
RTM is currently one of the best migration methods. It has excellent performance in complex structure imaging and amplitude preservation. In this section, we use the Marmousi model to compare our proposed method with the RTM method. In this experiment, we use a model with a flat surface because we want to eliminate the influence of terrain on the imaging effect and focus on the imaging differences between the two methods.
The Marmousi model is shown in Figure 4. The model space is of 1200 560 in grid size ( m). We place sources at shot locations and m. The source wavelet is a Ricker of 30 Hz central frequency. For each shot, the receivers are put at and m. The imaging results of RTM and our proposed method are shown in Figure 5a  From Figure 2, we can see that the imaging effects of the three methods are very close, and the position of the reflective layer is accurate. In Figure 3, we find that the calculated reflection coefficient curves of Method 2 and 3 are in good agreement with their corresponding theoretical curves, but the calculated reflection coefficient curve of Method 1 is inconsistent with the theoretical curve. Combining Figure 3 and Table 1, we can quantitatively analyze that the relative error of Method 1 is larger than that of the other two methods, and the performance of Method 2 and Method 3 are relatively close. These results indicate that Method 2 and 3 have better amplitude-preserving performance in this experiment. Table 1. The mean-square error of the relative errors (400-1400 m).

Method 1 Method 2 Method 3
Mean-square error 0.0384 0.0283 0.0274 In this experiment, there is a lateral velocity change in the upper layer of the velocity model. The vertical partial derivative calculation method used in Method 1 needs to assume a uniform surface velocity. When this assumption is not satisfied, there will be a non-negligible error in the calculated derivative. It will affect the amplitude-preserving performance of the imaging result. The performance of Method 2 and Method 3 are very close, and both methods can handle lateral surface velocity variation. However, these two methods are based on different acquisition systems: Method 2 is based on a dual-sensor acquisition system, and Method 3 is based on a conventional acquisition system.
All in all, compared with traditional methods, our proposed vertical partial derivative calculation method can handle lateral surface velocity variation. Compared with the method based on a dual-sensor acquisition system, our imaging effect is similar, but our acquisition system is simpler and more practical.

Marmousi Model
RTM is currently one of the best migration methods. It has excellent performance in complex structure imaging and amplitude preservation. In this section, we use the Marmousi model to compare our proposed method with the RTM method. In this experiment, we use a model with a flat surface because we want to eliminate the influence of terrain on the imaging effect and focus on the imaging differences between the two methods.
The Marmousi model is shown in Figure 4. The model space is of 1200 × 560 in grid size (dx = dz = 6.25 m). We place sources at shot locations x s i = i * 4dx, i = 0, 1, . . . , 239 and z = 0 m. The source wavelet is a Ricker of 30 Hz central frequency. For each shot, the receivers are put at x r j,i = x s i + jdx, j = 0, 1, . . . , 239 and z = 0 m. The imaging results of RTM and our proposed method are shown in Figure 5a,b, respectively. For a better comparison, we zoom in and display part of the imaging result in Figures 6 and 7. Figure 6 corresponds to the black rectangular part in Figure 4. Figure 7 corresponds to the white rectangular part. In order to compare the amplitude-preserving performance of the two methods, we extract the imaging amplitudes of the two logs (the location of log L1 and log L2 is shown in Figure 4) and compare them with the corresponding theoretical amplitudes. The theoretical amplitude here is obtained by convolution of reflection coefficient and wavelet. Figure 8 is the result of log L1. Figure 9 is the result of log L2.
Appl. Sci. 2021, 11, 3010 8 of 14 white rectangular part. In order to compare the amplitude-preserving performance of the two methods, we extract the imaging amplitudes of the two logs (the location of log L1 and log L2 is shown in Figure 4) and compare them with the corresponding theoretical amplitudes. The theoretical amplitude here is obtained by convolution of reflection coefficient and wavelet. Figure 8 is the result of log L1. Figure 9 is the result of log L2. From Figure 5, we can see that both methods can image the model well, but our proposed method is better than the RTM method in imaging of deep structures. Starting from a depth of about 3000 m, the imaging amplitude of the RTM method decreases significantly as the depth increases, and even some structures are missing (for example, the yellow rectangular area in Figure 5). For medium depth, we can see from Figures 6 and 7 that the imaging results of our method have better details than those of the RTM method. In Figures 8 and 9, we can find that the imaging amplitude curve of the RTM method agrees better with the theoretical curve than our method in the shallow part but at medium and deep depths, our method has obviously better performance than RTM. We can intuitively see that the imaging amplitude of the RTM will decrease as the depth increases, but our proposed method can maintain the amplitude very well. It proves that our method has better amplitude-preserving performance. This is why the imaging effect of our method is better than that of RTM in the deep part.                From Figure 5, we can see that both methods can image the model well, but our proposed method is better than the RTM method in imaging of deep structures. Starting from a depth of about 3000 m, the imaging amplitude of the RTM method decreases significantly as the depth increases, and even some structures are missing (for example, the yellow rectangular area in Figure 5). For medium depth, we can see from Figures 6 and 7 that the imaging results of our method have better details than those of the RTM method. In Figures 8 and 9, we can find that the imaging amplitude curve of the RTM method agrees better with the theoretical curve than our method in the shallow part but at medium and deep depths, our method has obviously better performance than RTM. We can intuitively see that the imaging amplitude of the RTM will decrease as the depth increases, but our proposed method can maintain the amplitude very well. It proves that our method has better amplitude-preserving performance. This is why the imaging effect of our method is better than that of RTM in the deep part.

Canadian Foothills Model
The Canadian Foothills model (shown in Figure 10) is a very classic SEG (Society of Exploration Geophysicists) model, which is often used to test the imaging performance of the migration algorithm in irregular topography scenario. This model shows the geological features from the Canadian Foothills in northeastern British Columbia [17]. The number of shots is 277, and the shot interval is 90 m. The first shot is placed at x = 0 m. There are 480 traces per shot. The offsets range is −3600 m to 3600 m for all shots except roll-in and roll-out. Figure 11 presents the topographic profile. The imaging result is shown in Figure 12. In this experiment, we select three logs (the location of log L3, log L4, and log L5 is shown in Figure 10). The comparison of the imaging amplitude and the theoretical amplitude of these three logs is shown in Figure 13.
From Figures 10 and 11, we can see that the surface conditions of this model are very complicated. The height difference between the highest point and the lowest point on the surface is 1537 m and there is a severe lateral velocity variation on the surface. In Figure 12, we can see that the imaging result can describe well the structures. All the faults are clearly visible. In Figure 13 we can find that the imaging amplitude curves are in good agreement with their theoretical curves. In this experiment, the deepest reflector is almost 9 km deep but our method can still accurately maintain the amplitude. These results prove that our method has good amplitude-preserving performance.

Canadian Foothills Model
The Canadian Foothills model (shown in Figure 10) is a very classic SEG (Society of Exploration Geophysicists) model, which is often used to test the imaging performance of the migration algorithm in irregular topography scenario. This model shows the geological features from the Canadian Foothills in northeastern British Columbia [17]. The number of shots is 277, and the shot interval is 90 m. The first shot is placed at m. There are 480 traces per shot. The offsets range is −3600 m to 3600 m for all shots except roll-in and roll-out. Figure 11 presents the topographic profile. The imaging result is shown in Figure 12. In this experiment, we select three logs (the location of log L3, log L4, and log L5 is shown in Figure 10). The comparison of the imaging amplitude and the theoretical amplitude of these three logs is shown in Figure 13.
From Figures 10 and 11, we can see that the surface conditions of this model are very complicated. The height difference between the highest point and the lowest point on the surface is 1537 m and there is a severe lateral velocity variation on the surface. In Figure 12, we can see that the imaging result can describe well the structures. All the faults are clearly visible. In Figure 13 we can find that the imaging amplitude curves are in good agreement with their theoretical curves. In this experiment, the deepest reflector is almost 9 km deep but our method can still accurately maintain the amplitude. These results prove that our method has good amplitude-preserving performance.

Discussion
Full-wave-equation depth extrapolation for migration is based on the two-way full-wave equation, therefore, the wave propagation can be described more accurately. However, because the conventional acquisition system cannot provide the vertical partial derivative of the surface wavefield, this method is very limited. The traditional method of obtaining the derivative requires the assumption that the surface velocity is constant and it is currently the most common method. This assumption can be satisfied

Discussion
Full-wave-equation depth extrapolation for migration is based on the two-way full-wave equation, therefore, the wave propagation can be described more accurately. However, because the conventional acquisition system cannot provide the vertical partial derivative of the surface wavefield, this method is very limited. The traditional method of obtaining the derivative requires the assumption that the surface velocity is constant and it is currently the most common method. This assumption can be satisfied

Discussion
Full-wave-equation depth extrapolation for migration is based on the two-way fullwave equation, therefore, the wave propagation can be described more accurately. However, because the conventional acquisition system cannot provide the vertical partial derivative of the surface wavefield, this method is very limited. The traditional method of obtaining the derivative requires the assumption that the surface velocity is constant and it is currently the most common method. This assumption can be satisfied at sea, but it is often unsatisfied on land. When the surface velocity has a lateral variation, there will be a non-negligible error in the calculated derivative. It will affect the amplitude-preserving performance of the migration result. This is proven in our horizontal layer model experiment. The vertical partial derivative calculation method we proposed overcomes the problem in that traditional methods need to assume uniform surface velocity although our method does have a disadvantage.For the reflected wave signal from the shallow part, the vertical partial derivative calculated by our method will have a certain error. From the results of the Marmousi model experiment, we can find that our method is slightly worse than the RTM method in terms of imaging amplitude accuracy in the shallow part (as shown in Figure 8, approximately at a depth of 500 m). However, as a whole, our method has obvious advantages in the description of structural details, the imaging of deep structures, and the ability to maintain amplitude. In actual seismic exploration, especially in petroleum seismic exploration, the depth of the target layer is generally several kilometers deep. Therefore, our method still has great practical application potential.

Conclusions
In this paper, we propose a new method for calculating the vertical partial derivative of the surface wavefield. It can be used in a topographically influenced scenario with lateral surface velocity variation. Our method solves the problem in that the traditional method needs to assume that the surface velocity is uniform. Based on the calculated derivative, we implement a depth-extrapolation-based full-wave-equation migration from topography by using direct downward continuation.
The results of the Marmousi model experiment show that our proposed method is superior to the conventional RTM algorithm in terms of imaging accuracy and amplitudepreserving performance at medium and deep depths. In the Canadian Foothills model experiment, we proved that our method can still accurately image complex structures and maintain amplitude in the irregular topography scenario with severe lateral velocity variation.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data used to support the findings of this study are available from the corresponding author upon request.