True Amplitude Angle Gathers from Reverse Time Migration by Wavefield Decomposition at Excitation Amplitude Time

Reservoir parameter estimation is one of the goals of amplitude-versus-angle (AVA) inversion and angle-domain common image gathers are the basis of AVA inversion. Therefore, the accuracy of kinematic and kinetic information on angle gathers is very important for reservoir characterization. Reverse time migration is one of the most physically accurate migration method. Generating angle gathers from reverse time migration with the Poynting vector method is very efficient. However, due to inaccurate angle measurement and uneven illumination, angle gathers calculated by the Poynting vector method are often not suitable for AVA inversion. In this paper, we propose an efficient method of angle gathers with accurate angular information and amplitude from reverse time migration. We firstly decompose source and receiver wavefield to their up-going and down-going parts by using analytic wavefield. We calculate propagation directions for source down-going wavefield and receiver up-going wavefield by the Poynting vector method and form the angle gathers with these angle information and decomposed wavefield. To reduce memory storage and improve computational efficiency, we decompose wavefield at excitation amplitude time by using a local spatial Fourier transform. We also use a spatial smoothed Poynting vector to improve the stability of angle measurement. We apply an illumination compensation image condition to recover the true amplitude. Numerical examples on Marmousi model and the SEAM two-dimensional (2D) model demonstrate the advantages of our proposed method. The angle gathers based on our method are cleaner with more focus on events energy and better continuity, suffering from less low-frequency noise in the shallow parts and with a distinct cutoff at large angle where reflection terminates. At last, we demonstrate the effectiveness of proposed method on a 2D marine field data example.


Introduction
The technology of amplitude versus angle (AVA) has been widely used to estimate subsurface petrophysical properties and predict fluid in the reservoirs for decades [1][2][3]. The input data for AVA analysis are angle gathers which are generated during the seismic migration process. Traditional approaches of calculating angle gathers are commonly implemented in ray-based migration methods since the ray-based propagator naturally provides angle information. However, the ray-based methods are limited by the high-frequency asymptotic approximation [4,5], thus, the wave equation based methods [6] were adopted for the angle gathers calculation.
As the wave equation-based methods do not provide propagation direction explicitly, several methods were proposed to extract angle information from wavefield. These methods can be summarized into two main categories. The first group applies the integral method, e.g., the local-plane wave decomposition and wavelet transform based methods [7][8][9], local Fourier transform based methods [10,11]. The second group applies differential methods, e.g., the Poynting vector or gradient vector-based methods [12][13][14][15][16]. Since the differential type methods only need to calculate spatial and temporal derivatives at specified locations, these methods are very efficient and easy to implement. Instead, the integral type methods need to analyze the local wavefield within a neighboring area, and these methods involve more computations, which are more time-consuming. However, compared to the differential methods, the integral methods are more stable, reliable and can handle the case of multi-arrivals.
To improve the accuracy and reliability of the differential methods, several methods were proposed, such as, stabilizing angle measurement in the time window [17], using time-shift Poynting vectors [18], smoothing the Poynting vectors in the space domain [13], using a least squares solution over the time or space [14], using the optical flow method to stabilize Poynting vector directions [19][20][21], using the traveltime gradient of first arrival to calculate the source-side wave propagation direction [22,23], or calculating angle just at most energetic arrival [16]. All these methods are effective to improve accuracy, and cannot deal with the situation of multiarrivals and overlapping events. Poynting vector methods were also applied to seismic illumination and resolution analyses, which can be used to correct the migration image [16]. The amplitude of angle gathers can be preserved by source illumination [24].
In this paper, we propose a new method to generate angle gathers under the reverse time migration frame, in which the angle information is calculated by a hybrid method. We first construct analytic wavefield and then decompose them into up-going and down-going parts in time-wavenumber domain. Then, we use the spatial weighted Poynting vectors method to calculate the wave propagation angle and obtain angle gathers with an angle domain illumination compensation imaging condition. Finally, we compare the angle gathers extracted by different methods and validate our method in the Marmousi model and 2D SEAM model to clarify the superiority of the proposed method.

Wavefield Decomposition in t-k Domain
Extracting the angle information by using the differential type methods (e.g., Poynting vectors) is highly efficient. However, these methods can only give one direction per spatial grid point, which means that if there are more than one wavefront coexists at one space point at the same time, the correct direction cannot be estimated with these methods. The overlapping of wavefield mainly occurs in two situations, the first situation is that the incident wave, reflected wave and transmitted wave coexist on the reflector, and the second one refers to that the multipath wavefield overlap at the same time. The former has great influence on the angle gathers, because it occurs on every reflector that needs to be imaged.
In order to deal with the overlapping wavefield at the reflector, we can separate source and receiver wavefield into their up-going, and down-going parts, respectively, where s(t, x), r(t, x) are source wavefield and receiver wavefield respectively. s u,d (t, x) are source up-and down-going parts respectively. r u,d (t, x) are receiver up-and down-going parts respectively. And x = (x, z), namely, s(t, x) = s(t, x, z) and r(t, x) = r(t, x, z). The conventional procedure to obtain the wave components in Equations (1) and (2) can be illustrated as follows [25]: where S(ω, k z ) is two-directional (2D) Fourier transform of s(t, z), ω is the angular frequency, k z is the vertical wavenumber, S U (ω, k z ) and S D (ω, k z ) are source up/down going wavefield in ω − k z domain respectively, likewise as receiver wavefield. The conventional method requires both signs of frequency and wavenumber to determine the up-going or down-going directions. Therefore, we need to store the entire wavefield, because the frequency information is obtained by performing Fourier transform on time axis, which is the slowest dimension of time-domain RTM. In order to reduce I/O cost, with the help of the characteristic that the analytic signal only has positive frequency components, we can construct analytic wavefield at first, and then decompose the analytic wavefield into up-going and down-going components through on-directional (1D) Fourier transform [26], whereŜ(t, k z ) is 1D Fourier transform of the analytic source wavefieldŜ(t, x), and k z is vertical wavenumber.Ŝ U,D (t, k z ) and s u,d (t, x) are source up/down-going components in time-wavenumber domain, and time-space domain, respectively.
operator takes out the real part of analytic wavefield. Likewise as receiver side wavefield. In order to avoid the convolution along the time axis of the Hilbert transform, we can apply the Hilbert transform to the source term and solve the wave equation to get the imaginary part of analytic wavefield [27]. However, decomposing wavefield by Fourier transform at each time step during wave extrapolation is very time-consuming especially in 3D case. Therefore, we just decompose wavefield in a local spatial window at excitation amplitude time [28,29] and separate wavefield into local angle components by weighted Poynting vectors [13]. The excitation time can be calculated by two methods: Ray tracing [30] and measuring the maximum amplitude of source wavefield [28,29]. We applied the second method, whereby the most energetic arrival time is used. We scan the source wavefield twice, first scan to locate the arrival time of most energetic amplitude, then decompose wavefield in the second scan. The twice-scanning procedure does not need additional computation since the source wavefield needs to be stored reproduced to correlate with receiver wavefield in RTM imaging.

Imaging in Angle Domain
Zero-lag cross correlation imaging condition can obtain RTM image stably. However, the resulting image does not have any physical meaning since the amplitude is unrelated to reflectivity. Deconvolution imaging condition can provide images related to reflectivity but becomes unstable when the denominator is small. In order to obtain stable angle gathers related to angle-dependent reflection coefficients, we apply a source-normalized imaging condition. The local angle image for an image location x can be obtained by, where s d (x, θ s , t) is the angle component of source-side down-going wavefield and r u (x, θ g , t) is the angle component of receiver-side up-going wavefield. θ s and θ g are local incident angle and local scattering angle respectively. We can convert the local angle image from (θ s ,θ g ) domain to (θ r ,θ d ) domain by a coordinate transform, θ d = (θ s + θ g )/2 (10) and θ r = (θ s − θ g )/2 (11) where θ r is the reflection angle and θ d is the target dipping angle. Since the AVA response is reflection-angle-dependent, the angle gathers can be obtained from a local angle image:

Angle Mesurement Analysis
We design a simple two layer model which is shown in Figure 1a to compare different angle measurement methods. We calculate angle gathers of one point (x = 3 km) with one shot (x = 2.4 km) by using different angle measurement methods. As mentioned above, the coexistence of incident wave, reflected wave and transmitted wave on the reflector will lead to the inaccurate angle measurement. In order to evaluate the influence of this situation on the angle measurement, we introduce a constant velocity model (Figure 1b) as migration velocity model which does not produce reflections. We use four methods to analyze the angle measurement, which are Poynting vector method, weighted Poynting vectors method, optical flow method and proposed method (wavefield decomposition plus weighted Poynting vectors). For simplicity, these four methods are referred to as Method 1, to Method 4, respectively. scattering angle respectively. We can convert the local angle image from ( domain by a coordinate transform, and where r  is the reflection angle and d  is the target dipping angle. Since the AVA response is reflection-angle-dependent, the angle gathers can be obtained from a local angle image:

Angle Mesurement Analysis
We design a simple two layer model which is shown in Figure 1a to compare different angle measurement methods. We calculate angle gathers of one point (x = 3 km) with one shot (x = 2.4 km) by using different angle measurement methods. As mentioned above, the coexistence of incident wave, reflected wave and transmitted wave on the reflector will lead to the inaccurate angle measurement. In order to evaluate the influence of this situation on the angle measurement, we introduce a constant velocity model ( Figure 1b) as migration velocity model which does not produce reflections. We use four methods to analyze the angle measurement, which are Poynting vector method, weighted Poynting vectors method, optical flow method and proposed method (wavefield decomposition plus weighted Poynting vectors). For simplicity, these four methods are referred to as Method 1, to Method 4, respectively. At first, we use the constant velocity model as migration velocity, and the angle gathers calculated by four methods are shown in Figure 2. Given there is a flat reflector at 1 km depth, the theoretical reflection angle is about 31° marked by the black solid line. As there is no interference of reflected wave, accurate reflection angle can be obtained with all these four methods. However, there are some smearing artifacts caused by unstable angle measurement in Figure 2a, which is the    Furthermore, we use true velocity model to generate the multi-shot angle gathers by using these four methods. The acquisition system is composed of 13 shots from 2.4 km to 3 km with a source interval of 50 m, and the angle gathers are extracted at 3 km, in which the theoretical reflection angle is about 0~31°. Figure 4a,b are the angle gathers calculated by Method 1 to 4. As illustrated in these figures, the angle gathers calculated by the first three methods are affected by smearing artifacts and low frequency noises (indicated by the red dashed box), while relatively accurate angle and sharp cutoff, at the end of reflection, can be obtained with Method 4. Moreover, smearing artifacts and low frequency noises can be efficiently eliminated by using Method 4.   Furthermore, we use true velocity model to generate the multi-shot angle gathers by using these four methods. The acquisition system is composed of 13 shots from 2.4 km to 3 km with a source interval of 50 m, and the angle gathers are extracted at 3 km, in which the theoretical reflection angle is about 0~31°. Figure 4a,b are the angle gathers calculated by Method 1 to 4. As illustrated in these figures, the angle gathers calculated by the first three methods are affected by smearing artifacts and low frequency noises (indicated by the red dashed box), while relatively accurate angle and sharp cutoff, at the end of reflection, can be obtained with Method 4. Moreover, smearing artifacts and low frequency noises can be efficiently eliminated by using Method 4. Furthermore, we use true velocity model to generate the multi-shot angle gathers by using these four methods. The acquisition system is composed of 13 shots from 2.4 km to 3 km with a source interval of 50 m, and the angle gathers are extracted at 3 km, in which the theoretical reflection angle is about 0~31 • . Figure 4a,b are the angle gathers calculated by Method 1 to 4. As illustrated in these figures, the angle gathers calculated by the first three methods are affected by smearing artifacts and low frequency noises (indicated by the red dashed box), while relatively accurate angle and sharp cutoff, at the end of reflection, can be obtained with Method 4. Moreover, smearing artifacts and low frequency noises can be efficiently eliminated by using Method 4.

Marmousi Model
In the first example, we extract angle gathers on the Marmousi model. The acquisition system is composed of 225 shots with a source interval of 50 m, and there are 200 receivers on the left side in each shot. The minimum and maximum offsets are 0 and 4975 m and the receiver interval is 25 m. Both the horizontal and vertical grid intervals are 12.5 m. The time interval is 1 ms. The source is a 20 Hz Ricker wavelet. In order to compare the effects of different methods, we first calculate the angle gathers, using four methods, based on the true velocity model. Partial results (indicated by the red dashed rectangle in Figure 5) are shown in Figure 6. Then, we use an incorrect velocity model as the migration velocity model to compute the angle gathers which are shown in Figure 7. As illustrated in Figure 6a-c, the angle gathers are very noisy and the events energy are not concentrated due to the inaccurate angle measurement. Methods 1 to 3 cannot handle with back-scattering artifacts if there are complex structures. On the contrary, the angle gathers in Figure 6d are with higher signal-to-noise ratio and clear cutoff at large angle which is shown in the red solid rectangles. Moreover, the shallow reflector marked in red dashed box can be clearly imaged in Figure 6d, while the obvious events cannot be obtained by the first three methods, especially at small angles.
In Figure 7, due to the incorrect migration velocity model, most of the source wavefields and receiver wavefields are not in the same wavepath. Therefore, there are fewer back-scattering artifacts in Figure 7a-c than in Figure 6a-c. The angle gathers in Figure 7d are with the best quality among the four results. In details, the angle gathers in Figure 7d are cleaner, more concentrated, with better continuity and sharper cutoff. It is also shown that the angle gathers calculated by Method 4 are more suitable for the migration velocity analysis because the curvatures are easier to pick up.
In Figure 8, we generate the stack images with the reflection angle from 0 to 50 • of four methods. Since large reflection angles are excluded, all images do not suffer from low-frequency noise. Compared with Figure 8a-c, the stack image by Method 4 in Figure 8d is with consistent layered structures in the shallow part marked by red dash rectangles, more balanced amplitudes, and sharper images in general. The areas marked by black square in Figure 8a,d are chosen to quantitatively compare the resolution of stack images. Figure 8e,f are the local image and wavenumber spectrum of stack image by Method 1, Figure 8g,h are the local image and wavenumber spectrum of stack image by Method 4. In Figure 8f, the energy is mainly concentrated on the low wavenumber part since the reflection angle is incorrectly estimated. After accurate angle calculation by Method 4, the energy distribution on the wavenumber spectrum is more balanced which can be shown in Figure 8h.

Marmousi Model
In the first example, we extract angle gathers on the Marmousi model. Hz Ricker wavelet. In order to compare the effects of different methods, we first calculate the angle gathers, using four methods, based on the true velocity model. Partial results (indicated by the red dashed rectangle in Figure 5) are shown in Figure 6. Then, we use an incorrect velocity model as the migration velocity model to compute the angle gathers which are shown in Figure 7. As illustrated in Figure 6a-c, the angle gathers are very noisy and the events energy are not concentrated due to the inaccurate angle measurement. Methods 1 to 3 cannot handle with back-scattering artifacts if there are complex structures. On the contrary, the angle gathers in Figure 6d are with higher signal-tonoise ratio and clear cutoff at large angle which is shown in the red solid rectangles. Moreover, the shallow reflector marked in red dashed box can be clearly imaged in Figure 6d, while the obvious events cannot be obtained by the first three methods, especially at small angles.    In Figure 7, due to the incorrect migration velocity model, most of the source wavefields and receiver wavefields are not in the same wavepath. Therefore, there are fewer back-scattering artifacts in Figure 7a-c than in Figure 6a-c. The angle gathers in Figure 7d are with the best quality among the four results. In details, the angle gathers in Figure 7d are cleaner, more concentrated, with better continuity and sharper cutoff. It is also shown that the angle gathers calculated by Method 4 are more suitable for the migration velocity analysis because the curvatures are easier to pick up.
In Figure 8, we generate the stack images with the reflection angle from 0 to 50° of four methods. Since large reflection angles are excluded, all images do not suffer from low-frequency noise. compare the resolution of stack images. Figure 8e,f are the local image and wavenumber spectrum of stack image by Method 1, Figure 8g,h are the local image and wavenumber spectrum of stack image by Method 4. In Figure 8f, the energy is mainly concentrated on the low wavenumber part since the reflection angle is incorrectly estimated. After accurate angle calculation by Method 4, the energy distribution on the wavenumber spectrum is more balanced which can be shown in Figure 8h.

SEAM 2D Model
In order to evaluate the quality of subsalt imaging, we investigated a 2D section of the SEAM phase I model in this example. The P wave velocity model is shown in Figure 9. The acquisition system is composed of 400 shots with a source interval of 80 m. The minimum and maximum offsets are 0 and 7980 m and the receiver interval is 40 m. The horizontal and vertical grid intervals are 20 m, and 10 m, respectively. The time interval is 1 ms. The source is a 20 Hz Ricker wavelet.

SEAM 2D Model
In order to evaluate the quality of subsalt imaging, we investigated a 2D section of the SEAM phase I model in this example. The P wave velocity model is shown in Figure 9. The acquisition system is composed of 400 shots with a source interval of 80 m. The minimum and maximum offsets are 0 and 7980 m and the receiver interval is 40 m. The horizontal and vertical grid intervals are 20 m, and 10 m, respectively. The time interval is 1 ms. The source is a 20 Hz Ricker wavelet. We selected a subsalt area marked by red dashed rectangles in Figure 9 to demonstrate the angle gathers. The results are shown in Figure 10. The angle gathers that were calculated by Method 2 (Figure 10b) and Method 3 ( Figure 10c) have a sharper cutoff at the large reflection angle than the angle gathers calculated by Method 1 (Figure 10a) because weighted Poynting vector and optical flow can stabilize the angle measurement. Nevertheless, these three methods do not separate the up-going and down-going wavefield, so the artifacts and low-frequency noises caused by backscattering wave cannot be effectively eliminated, and the errors of angle measurement still exist. By contrast, Method 4 decomposes wavefield into up-going and down-going parts, and is effective in eliminating the influence of multiarrivals and obtains more accurate angle measurement. As a result, angle gathers in Figure 10d are with higher signal-to-noise ratio, more balanced amplitude, and wider angle coverage. We selected a subsalt area marked by red dashed rectangles in Figure 9 to demonstrate the angle gathers. The results are shown in Figure 10. The angle gathers that were calculated by Method 2 (Figure 10b) and Method 3 ( Figure 10c) have a sharper cutoff at the large reflection angle than the angle gathers calculated by Method 1 (Figure 10a) because weighted Poynting vector and optical flow can stabilize the angle measurement. Nevertheless, these three methods do not separate the up-going and down-going wavefield, so the artifacts and low-frequency noises caused by backscattering wave cannot be effectively eliminated, and the errors of angle measurement still exist. By contrast, Method 4 decomposes wavefield into up-going and down-going parts, and is effective in eliminating the influence of multiarrivals and obtains more accurate angle measurement. As a result, angle gathers in Figure 10d are with higher signal-to-noise ratio, more balanced amplitude, and wider angle coverage.

Field Data
Finally, we test our method with a 2D marine field data. Figure 11 shows the migration velocity model. The acquisition system is composed of 530 shots with a source interval of 50 m, and there are

Field Data
Finally, we test our method with a 2D marine field data. Figure 11 shows the migration velocity model. The acquisition system is composed of 530 shots with a source interval of 50 m, and there are 360 receivers on the left side in each shot and the receiver interval is 12.5 m. The minimum and maximum offsets are 183 and 4670.5 m. We generated the angles using Method 1 and Method 4 shown in Figure 12a,b. We can see that the angle gathers calculated by Method 4 are more energy-balanced and continuous. We also generated the stack images with the reflection angle from 0 to 50 • by Method 1 and Method 4 shown in Figure 13a,b respectively. The stack image by Method 4 is of more balanced amplitude, and more sharper image. The low wavenumber artifacts in shallow parts of Figure 13a are also eliminated in Figure 13b, because the proposed method can estimate the reflection angle more accurately.
Energies 2020, 13, x 12 of 16 maximum offsets are 183 and 4670.5 m. We generated the angles using Method 1 and Method 4 shown in Figure 12a,b. We can see that the angle gathers calculated by Method 4 are more energybalanced and continuous. We also generated the stack images with the reflection angle from 0 to 50° by Method 1 and Method 4 shown in Figure 13a,b respectively. The stack image by Method 4 is of more balanced amplitude, and more sharper image. The low wavenumber artifacts in shallow parts of Figure 13a are also eliminated in Figure 13b, because the proposed method can estimate the reflection angle more accurately. Figure 11. The migration velocity model from a marine field data.  Energies 2020, 13, x 12 of 16 maximum offsets are 183 and 4670.5 m. We generated the angles using Method 1 and Method 4 shown in Figure 12a,b. We can see that the angle gathers calculated by Method 4 are more energybalanced and continuous. We also generated the stack images with the reflection angle from 0 to 50° by Method 1 and Method 4 shown in Figure 13a,b respectively. The stack image by Method 4 is of more balanced amplitude, and more sharper image. The low wavenumber artifacts in shallow parts of Figure 13a are also eliminated in Figure 13b, because the proposed method can estimate the reflection angle more accurately. Figure 11. The migration velocity model from a marine field data.

Computational Issues
From the formulation of wave decomposition and angle domain imaging, we can see that intensive computations are required to decompose wavefield into up/dwon-going parts and calculate angle components of wavefield. By constructing the analytic wavefield and applying the Hilbert transform on source term, 2D Fourier transform can be simplified to 1D Fourier transform and convolution along time axis can be avoided in the process of wavefield decomposition. However, a large number of Fourier transform is still required in wavefield decomposition. In terms of the 2D model of N N  times 1D Fourier transform is required at each time during the wavefield extrapolation, and the number of samples for each Fourier transform is z N . Except for the excitation time, most of the wavefield are with minimal contribution to imaging. Therefore, we propose the scheme in which we only decompose the wavefield by using a local 1D Fourier transform at the excitation time, and the number of samples for each Fourier transform is 32, much less than z N . The excitation time for each spatial grid can be easily obtained by extrapolating the wavefield twice. In RTM imaging, the source wavefield needs to be either stored or regenerated to correlate with receiver wavefield. Therefore, twice-extrapolating procedure at source-side does not require additional computation.

Computational Issues
From the formulation of wave decomposition and angle domain imaging, we can see that intensive computations are required to decompose wavefield into up/dwon-going parts and calculate angle components of wavefield. By constructing the analytic wavefield and applying the Hilbert transform on source term, 2D Fourier transform can be simplified to 1D Fourier transform and convolution along time axis can be avoided in the process of wavefield decomposition. However, a large number of Fourier transform is still required in wavefield decomposition. In terms of the 2D model of N x × N z grids for example, the number of time sampling points during wavefield extrapolation is N t . Decomposing the wavefield into up/down-going parts according to the scheme proposed by Wang et al. (2016) N t × N x times 1D Fourier transform is required at each time during the wavefield extrapolation, and the number of samples for each Fourier transform is N z . Except for the excitation time, most of the wavefield are with minimal contribution to imaging. Therefore, we propose the scheme in which we only decompose the wavefield by using a local 1D Fourier transform at the excitation time, and the number of samples for each Fourier transform is 32, much less than N z . The excitation time for each spatial grid can be easily obtained by extrapolating the wavefield twice. In RTM imaging, the source wavefield needs to be either stored or regenerated to correlate with receiver wavefield. Therefore, twice-extrapolating procedure at source-side does not require additional computation.

Discussion
In the proposed method, we adopted the Poynting vector to measure the wave propagation direction. Compared to integral methods, the Poynting vector method was much more efficient, but less accurate. To avoid the situation that the incident wave and reflected wave coexist at the Energies 2020, 13, 6204 14 of 16 reflector, which will cause the incorrect angle measurement, we decompose the wavefield into its up-going and down-going parts. Since the wavefield of excitation time has the greatest contribution to imaging, we only decomposed the wavefield at the excitation time, which can also improve efficiency. To stabilize the angle measurement, a spatial weighted Poynting vector was conducted. An angle domain source illumination normalized imaging condition was applied to obtain the angle gathers. The numerical experiments on the Marmousi model demonstrate that the proposed method can improve the accuracy of angle measurement with the decomposition of up-going and down-going wavefield. When the velocity model is inaccurate, the curvatures of the angle gathers, calculated by the proposed method, are easier to pick up, indicating that the proposed method is possible to be applied in migration velocity analysis. The results of the SEAM 2D model show that the amplitude of subsalt areas can be effectively recovered with the proposed method through angle domain illumination compensation. The example on the field data also demonstrate that the proposed method can gives better quality angle gathers than the Poynting vector method. However, since only the primary wave is considered, and only two parts of up-going and down-going are decomposed, the effectiveness of this method, under a more complex environment, may need further investigation.

Conclusions
In this paper, an efficient method is proposed for extracting angle gathers from reverse time migration. At first, the wavefield is decomposed into its up-going and down-going components by constructing the analytic wavefield. To improve the efficiency, the wavefield is decomposed only at excitation time, thus, only the primary wave is contributed to imaging. Then, a weighted Poynting vector is adopt to calculate the wave propagation angle. A simple two-dimensional model was used to verify that the accuracy of angle measurement can be effectively improved by the proposed method. To obtain the angle gathers related to the angle-dependent reflection coefficients, an angle domain illumination normalized imaging condition was used as well. In this study, the present two numerical examples and a 2D marine field were used as examples to validate the effectiveness of the proposed method.