Radon Transform Based on Waveform for AVO-Preserving Data Construction

: The traditional hyperbolic Radon transform suffers from the major problem of how to both obtain a high resolution and preserve the amplitude variation with offset (AVO). In the Radon domain, high resolution (sparseness) is a valid criterion. However, if a sparse model is obtained in the Radon domain due to averaging along the offset direction, then it is not possible to preserve the AVO in the inversion data. In addition, hyperbolic Radon transform has a time-variant kernel based on a traditional iterative algorithm, the conjugate gradient (CG), which requires signiﬁcant computation time. To solve these problems, we propose a Radon transform based on waveform that contains both cycle and amplitude characteristics of seismic waves. The new transform entails creating an upper envelope for the seismic data and computing a preliminary forward Radon transform in the time domain. The forward Radon transform incorporates a priori information by measuring the energy of each slowness (p) trace to obtain the high-resolution result of the Radon domain. For AVO preserving, the proposed method uses polynomials to describe the AVO characteristics in the inverse Radon transform based on the least-squares inversion. Besides amplitude preserving and high resolution, the proposed method avoids using CG and greatly reduces the cost of computing hyperbolic Radon transform in the time domain. In applications to both synthetic and ﬁeld data, waveform Radon transform (WRT) has a better performance than the conjugate gradient Radon transform (CGRT).


Introduction
Over the past two decades, Radon transform has been widely used in seismic data processing, especially in multiple attenuation and seismic data interpolation [1][2][3][4][5][6][7]. Radon transform can be divided into three different types based on its shape: line, parabolic, or hyperbolic Radon transform. Thorson et al. [8] proposed a time-domain hyperbolic least-squares method, which led to the high-resolution transform at the expense of a very large matrix computation. Hampson [9] (1986) suggested a frequency domain least-squares parabolic transform, which is more efficient. Therefore, the frequency domain quickly became the industry standard. Sacchi and Ulrych [10] (1995) implemented the sparse Radon transform in the frequency domain, which obtained a high-resolution result. Their frequency domain technique is now generally used in seismic processing. Sparse Radon transform improves the resolution but introduces new problems that make it difficult to preserve the AVO.
When the resolution is high enough, there is only one sample in the Radon model associated with every reflection event, and the inverse transform amplitude is constant. Thus, the resolution and AVO characteristics are conflicting. Kabir [11] (1999) pointed out that Radon transform is not an orthogonal transform, and the inverse transform cannot completely reconstruct the original seismic data, especially the amplitude distortion at the near offset, resulting in the AVO analysis error. Johansen [12] (1995) used the orthogonal polynomials to fit the AVO characteristics of the CDP gathers, which showed that AVO can be characterized by fewer coefficients of low-order orthogonal polynomials. Mark, N [13] (2004) proposed a high-resolution Radon transform in the t − x domain using the "intelligent" prioritization of the Gauss-Seidel estimation sequence. Xue [14] (2009) proposed a de-multiple method based on Radon and orthogonal transform which showed a good performance in retaining the AVO properties. Wang [15] (2014) split the seismic gather into two Radon gathers, which reduced the amplitude distortion and preserved the AVO phenomena compared with the classical sparse Radon transform.
Another problem of high-resolution Radon transform is the large computation time. Compared with the frequency-domain Radon transform, the time-domain Radon transform achieves a higher resolution because the Radon model has a higher sparseness in the time domain. The downside of the time-domain Radon transform is its significantly increased computation burden. Trad et al. [16] (2003) proposed a method to compute the timeinvariant Radon transform in the time domain applying frequency domain operators. This approach combines the high resolution of the time-domain Radon transform and the computational efficiency of the frequency-domain Radon transform. However, this technique cannot implement the hyperbolic Radon Transform, which is time-variant. In general, a true hyperbolic Radon transform produces better results for data with large moveout than the parabolic version. In addition, the hyperbolic Radon transform can be directly applied to CMP data without NMO-corrected data. The large computation cost is the main problem of the hyperbolic Radon transform. Due to its time-variant kernel, it is not possible to use fast solutions instead of iterative methods such as the conjugate gradient, which increase the computation time [17][18][19][20][21][22][23].
This paper describes a technique for implementing a high-resolution hyperbolic Radon transform-preserving AVO. To obtain the high-resolution hyperbolic Radon transform, we fully consider the characteristics of the seismic data. Previous research on the Radon transform has mainly focused on the kinematic characteristics of the seismic data signal, rarely considering the waveform properties containing energy and cycle. This paper uses a new method to improve the resolution based on the waveform energy and cycle. First, the technique entails making an upper envelope for the seismic data and computing a preliminary forward Radon transform in the time domain. Next, the second-forward Radon transform incorporates a priori information by computing the energy of each p trace in the preliminary transform. We can obtain the high-resolution Radon model by twice-forward Radon transform. To address the problem of lacking AVO, we modify the Radon function by incorporating the least-squares principle in the inversion Radon transform using polynomials to describe the AVO characteristics. The proposed method avoids calculating the time-variant kernel. Compared with CGRT, the proposed method has a high computation efficiency.

Hyperbolic Radon Transform
In the layered subsurface model, the reflection events fit the hyperbolic trajectory in the common midpoint (CMP) gather, d(t, x). Therefore, the hyperbolic Radon transform can be used to process the CMP gather data. The transform formula is as follows where t is the zero-offset two-way time, x is the offset, m(τ, p) is the hyperbolic Radon domain, τ is intercept time and p is the parameter that defines the curvature. k = 1, 2 . . . K, K is the number of the offset x, n = 1, 2 . . . N, N is the number of the curvature p and i is the discretize sequence number for the time t. Generally, in seismic data processing, the hyperbolic Radon transform first defines the inverse transform and then obtains the forward transform formula using the least-squares principle. In matrix-vector notation, the inverse hyperbolic Radon transform formula can be written as: where d is the CMP gather, L is the forward hyperbolic Radon transform operator and m is the hyperbolic Radon model. The objective function of the least-squares solution is as follows: The least-squares solution of the forward hyperbolic Radon transform is: L T is the adjoint hyperbolic Radon operator. A sparse constraint in the time domain can improve the resolution of the Radon panel. To obtain a high-resolution result, the mixed normed is introduced to constrain the inversion. The optimization problem can be written as: where the l 2 norm constrains the misfit of the inversion, the l 1 norm constrains the sparsity of the solution and λ is the tradeoff parameter, which controls the sparsity of the Radon panel. If λ increases, then the sparsity of the Radon panel would be improved, that is, a high resolution of the Radon panel would be obtained. However, due to averaging along the offset direction, the sparsity of the Radon panel does not preserve the AVO.

Conjugate Gradient Method
The conjugate gradient method (CG) is usually used to solve Equation (6). The solving process can be understood as the problem of solving linear equations A·x = b using CG given the least-squares solution for the overdetermined problem and the minimum norm solution for the underdetermined problem. The steps of CG are given as follows, assuming four vector sequences r k , r k , p k , p k , k = 1,2..., and p 1 = r 1 , p 1 = r 1 : The large amount of computation is an inherent problem in solving the Radon transform in the time domain. The operator L and L T is the large sparse matrix. Each iteration requires the multiplication of the large sparse matrices L and L T with vectors, which requires significant calculation.

Waveform Radon Transform
We propose a new method based on the waveform to realize the hyperbolic Radon transform. First, we provide an interpretation of low resolution in Radon transform in the time domain. We simply illustrate the forward transform process in Figure 1. In Figure 1,

Waveform Radon Transform
We propose a new method based on the waveform to realize the hyperbolic Radon transform. First, we provide an interpretation of low resolution in Radon transform in the time domain. We simply illustrate the forward transform process in Figure 1. In Figure 1,  To obtain a high resolution, it is necessary to reduce the influence of the wrong stacking trajectory. Therefore, in the forward transform, the ideal situation is for all stacking trajectory to match the reflected waves. To achieve this, we introduced the information of the energy point of each curvature in the Radon panel and the width of the wave cycle. The main difference between the correct trajectory and the wrong trajectory is their corresponding energy points in the Radon domain: the correct trajectory's energy point is larger, and the shape is closer to the waveform. Besides, the curvature of each cycle of the waveform is consistent, with every cycle corresponding a correct curvature.
However, some factors that influence the accuracy of the energy point, the noise and the negative values of the wavelet, such as the trough. In addition, the waveform has a complex shape, so the wave cycle can be difficult to obtain. To overcome these problems, we used the waveform mode decomposition (WMD) to process the wavelet before the Radon transform [24]. We used the upper envelope in the WMD to present the wavelet. Taking the Rick wavelet as an example, the waveform and the upper envelope are shown in Figure 2. The upper envelope can reflect the position of the wave crest and the wave cycle, . To obtain a high resolution, it is necessary to reduce the influence of the wrong stacking trajectory. Therefore, in the forward transform, the ideal situation is for all stacking trajectory to match the reflected waves. To achieve this, we introduced the information of the energy point of each curvature in the Radon panel and the width of the wave cycle. The main difference between the correct trajectory and the wrong trajectory is their corresponding energy points in the Radon domain: the correct trajectory's energy point is larger, and the shape is closer to the waveform. Besides, the curvature of each cycle of the waveform is consistent, with every cycle corresponding a correct curvature.
However, some factors that influence the accuracy of the energy point, the noise and the negative values of the wavelet, such as the trough. In addition, the waveform has a complex shape, so the wave cycle can be difficult to obtain. To overcome these problems, we used the waveform mode decomposition (WMD) to process the wavelet before the Radon transform [24]. We used the upper envelope in the WMD to present the wavelet. Taking the Rick wavelet as an example, the waveform and the upper envelope are shown in Figure 2. The upper envelope can reflect the position of the wave crest and the wave cycle, t w . Therefore, obtain high-resolution results in the Radon domain, two Radon transformations are needed. In the first Radon transform, the correct stacking trajectory are obtained using the energy points in the Radon domain. In the second transform, the − data are stacked along the correct stacking trajectory selected in the first transform. The specific steps are as follows: 1. Make the first forward hyperbolic Radon transform in the time domain: stack the cur- Arrange the samples in the Radon panel in order of large energy to small energy and obtain an array S = S 1 , S 2 , . . . . . . S nτ×np . Each sample in array S has an intercept time τ and curvature p. The sample S 1 is the largest energy sample in the Radon model and corresponds to τ S 1 , p S 1 . S n (n = 1, 2 . . . . . . nτ × np) corresponds to (τ S n , p S n ).

3.
Make the second forward Radon transform in the time domain and the transform sequence accord to array S. First, transform the τ S 1 , p S 1 in the t − x data gather and eliminate the amplitude energy of the stacking path. Use p S 1 to transform and eliminate these trajectories. Then, proceed in sequence by array S. The specific formula is as follows: n = 1, 2 . . . . . . nτ × np, K is the number of the offset x and d S n (t, x) is the t − x data gather when transform the S n .

4.
Loop nτ × np times in the order of array S. Every loop eliminates a reflection event trace. After loops are completed, we obtain the high-resolution Radon model. Usually, it does not need to Loop nτ × np times, as in steps 2-3. Finishing nτ loops is enough. Because of the number of the reflection event traces, there are no more than nτ sampling points.
After completing the twice-forward Radon transform as described above, we obtained a high-resolution solution in the Radon panel. However, the proposed hyperbolic Radon transform based on the principle of energy ordering can obtain the spare solution. The AVO of the original seismic data cannot be preserved using the traditional inverse Radon transform.

Inversion Transform Preserving AVO
To preserve the AVO, we incorporated the least-squares principle when designing the inverse Radon transform. Unlike the traditional Radon transform, the inverse Radon transform defines the inverse transform first and uses the least-squares principle to obtain the forward transform formula. The proposed method already obtained the high-resolution Radon model result based on the energy principle, so we used the least-squares result based on the energy principle, to obtain the inverse transform. The misfit between the original seismic data and inverse seismic data was minimized and subjected to the least-squares principle.
The AVO characteristics of the CDP gathers can be fitted by polynomials [25]. In addition, fewer coefficients of the low-order polynomials can describe the main energy of the amplitude (Johansen et al., 1995 [12]). Therefore, the AVO characteristics of the seismic reflection wave can be described as: F(x) is the amplitude energy, x is the offset, and c 0 , c 1 and c 2 are the coefficients. The formula of the inverse transform is as follows: is the inverse Radon transform result of proposed method. The objective function of the least square is as follows: By solving Equation (13), we obtain the values of c 0 , c 1 and c 2 .
The whole process of the proposed method does not include the time-variant kernel. It clearly improves the computational efficiency.

Synthetic Examples
To evaluate the performance of the proposed method, we designed synthetic data with a Ricker wavelet which showed strong amplitude variation with offset. Figure 3a shows the synthetic CMP gather containing 70 offsets from 0 to 2760 m. There were three reflecting waves located at 500 ms, 900 ms and 1200 ms, respectively. Two the reflecting waves were variable amplitude events at 500 ms and 900 ms, respectively. Figure 3b shows the Ricker wavelet and its upper envelope of the first trace and the last trace in Figure 1a. To accurately describe the AVO characteristics the AVO characteristics of the reflections located at 500 ms and 900 ms are provided in Figure 4.    To test the data reconstruction performance, we created some empty channels in the synthetic data, as shown in Figure 5a. We used the waveform Radon transform (WRT) and CGRT to reconstruct the data gather. The Radon panel results in Figure 5b,c show that the resolution of the WRT in the Radon panel was higher than the CGRT, and the WRT waveform was closer to the Ricker wavelet. Figure 6 shows the inverse transform result and its residual.  To test the data reconstruction performance, we created some empty channels in the synthetic data, as shown in Figure 5a. We used the waveform Radon transform (WRT) and CGRT to reconstruct the data gather. The Radon panel results in Figure 5b,c show that the resolution of the WRT in the Radon panel was higher than the CGRT, and the WRT waveform was closer to the Ricker wavelet. Figure 6 shows the inverse transform result and its residual.  Compared with the CGRT, the WRT inverse transform results recovered the original t − x data better, as shown in Figure 6. Figure 7 shows the AVO characteristic between the true and predicted data using the CGRT and proposed WRT. It shows that the WRT was able to preserve the AVO on the full range of offsets. Figure 8 shows the waveform of the original data at the first trace and one of the reconstruction traces. The black line represents the original data, the red line and the blue line represent the reconstruction data of the CGRT and WRT, respectively. Compared with the CGRT, the WRT kept an accurate waveform with little error. The waveform error of the WRT increased at the reconstruction trace. However, the WRT still kept an accurate waveform compared with the CGRT.  Compared with the CGRT, the WRT inverse transform results recovered the original − data better, as shown in Figure 6. Figure 7 shows the AVO characteristic between the true and predicted data using the CGRT and proposed WRT. It shows that the WRT was able to preserve the AVO on the full range of offsets.    Figure 8 shows the waveform of the original data at the first trace and one of the reconstruction traces. The black line represents the original data, the red line and the blue line represent the reconstruction data of the CGRT and WRT, respectively. Compared with the CGRT, the WRT kept an accurate waveform with little error. The waveform error of the WRT increased at the reconstruction trace. However, the WRT still kept an accurate waveform compared with the CGRT. In order to quantitatively describe the results of the data reconstruction, we use signal-to-noise ratio, defined as In order to quantitatively describe the results of the data reconstruction, we use signal-to-noise ratio, defined as SNR = 10 log 10 ||x|| 2 2 ||x −x|| 2 2 (19) where x represents original full data andx represents reconstructed data. The SNR of the reconstructed results obtained by the CGRT and the WRT are 10.49 dB and 22.69 dB, respectively. We recorded the running time of the CGRT and WRT. The hardware of the computer used in the experiment included an intel(R) Core (TM) i7-4790 cpu@3.60 GHz. The computational time of the CGRT was 135 s, and that of the WRT was 15 s. These results prove that the WRT was faster than the CGRT.

Field Data Examples
We used 2D field data to demonstrate the effectiveness of the WRT, as shown in Figure 9a. There were 126 offsets in total, the offset spacing was 24 m and each offset had 650 sampling points. The sampling interval was 4 ms, and the recording length was 2.6 s. Figure 9b shows the waveform and its upper envelope for trace 1 and trace 55.
To test the reconstruction performance, we designed empty regions (from trace 61 to trace 70) in the seismic data, as shown in Figure 10a. From Figure 10b,c the WRT was significantly better than the CGRT with respect to the resolution in the Radon panel. Then, we used the CGRT and WRT to reconstruct the field seismic data gather. The reconstruction result is shown in Figure 11, and the area between the two red lines is the reconstruction region. It is clear that the reconstruction region of the WRT was better than the CGRT. that the WRT was faster than the CGRT.

Field Data Examples
We used 2D field data to demonstrate the effectiveness of the WRT, as shown in Figure 9a. There were 126 offsets in total, the offset spacing was 24 m and each offset had 650 sampling points. The sampling interval was 4 ms, and the recording length was 2.6 s. Figure 9b shows the waveform and its upper envelope for trace 1 and trace 55. To test the reconstruction performance, we designed empty regions (from trace 61 to trace 70) in the seismic data, as shown in Figure 10a. From Figure 10b,c the WRT was significantly better than the CGRT with respect to the resolution in the Radon panel. Then, we used the CGRT and WRT to reconstruct the field seismic data gather. The reconstruction result is shown in Figure 11, and the area between the two red lines is the reconstruction region. It is clear that the reconstruction region of the WRT was better than the CGRT.   Figure 12 shows the waveform of the original data and the inversion results at trace 10 and trace 64 at the reconstruction region. The black line represents the original data, while the red line and the blue line represent CGRT and WRT, respectively. Compared with the CGRT, for the field data, the WRT kept the waveform more accurate and preserved the AVO. The SNR of the reconstructed results for the CGRT and for the WRT are 4.32 dB and 9.71 dB, respectively.   Figure 12 shows the waveform of the original data and the inversion results at tr 10 and trace 64 at the reconstruction region. The black line represents the original da while the red line and the blue line represent CGRT and WRT, respectively. Compar with the CGRT, for the field data, the WRT kept the waveform more accurate and p served the AVO. The SNR of the reconstructed results for the CGRT and for the WRT 4.32 dB and 9.71 dB, respectively. For the field data, the differences in the calculation time were more apparent. The running time of the CGRT was 2850.137 s, while that of the WRT was 114.257 s. This proves that WRT is at least one order of magnitude faster than CGRT.

Conclusions
In this paper, we proposed a waveform Radon transform (WRT) method considering both wave shape characteristics and energy characteristics of seismic data, which can preserve AVO. The proposed method used an upper envelope instead of the seismic waveform to obtain the energy and wave cycle information. The forward Radon transform incorporated the prior information of the waveform to improve the resolution in the Radon For the field data, the differences in the calculation time were more apparent. The running time of the CGRT was 2850.137 s, while that of the WRT was 114.257 s. This proves that WRT is at least one order of magnitude faster than CGRT.

Conclusions
In this paper, we proposed a waveform Radon transform (WRT) method considering both wave shape characteristics and energy characteristics of seismic data, which can preserve AVO. The proposed method used an upper envelope instead of the seismic waveform to obtain the energy and wave cycle information. The forward Radon transform incorporated the prior information of the waveform to improve the resolution in the Radon model. The inverse Radon transform based on the least-squares inversion was implemented in the time domain, incorporating the coefficients of low-order polynomials to describe the main energy of the amplitude, which can preserve a good AVO.
We compared the performance of traditional sparse Radon transform (CGRT) and wave Radon transform (WRT). WRT obtained a high resolution in the Radon panel, and the inverse transform was more similar to the original seismic data, preserving a good AVO. In addition, computational efficiency is an important advantage of WRT, especially for field data. In the future, WRT can be further applied for multiple scenarios.