TIME-INVARIANT RADON TRANSFORM BY GENERALIZED FOURIER SLICE THEOREM

. Time-invariant Radon transforms play an important role in many ﬁelds of imaging sciences, whereby a function is transformed linearly by inte- grating it along speciﬁc paths, e.g. straight lines, parabolas, etc. In the case of linear Radon transform, the Fourier slice theorem establishes a simple ana- lytic relationship between the 2-D Fourier representation of the function and the 1-D Fourier representation of its Radon transform. However, the theorem can not be utilized for computing the Radon integral along paths other than straight lines. We generalize the Fourier slice theorem to make it applicable to general time-invariant Radon transforms. Speciﬁcally, we derive an analytic expression that connects the 1-D Fourier coeﬃcients of the function to the 2-D Fourier coeﬃcients of its general Radon transform. For discrete data, the model coeﬃcients are deﬁned over the data coeﬃcients on non-Cartesian points. It is shown numerically that a simple linear interpolation provide satisfactory results and in this case implementations of both the inverse operator and its adjoint are fast in the sense that they run in O ( N log N ) ﬂops, where N is the maximum number of samples in the data space or model space. These two canonical operators are utilized for eﬃcient implementation of the sparse Radon transform via the split Bregman iterative method. We provide numer- ical examples showing high-performance of this method for noise attenuation and waveﬁeld separation in seismic data.

1. Introduction. The Radon transform has a long history in image processing and seismic data processing. In image processing it is often called the Hough transform [13] and it has been used for pattern recognition and feature extraction from images. The Radon transform has also been used in seismic data processing to solve wave propagation and inversion problems in media where the subsurface speed only varies with depth [28]. In the context of our work, we will discuss the Radon transform and a novel implementation of it, for seismic data processing. Radon transforms have been extensively used for random and coherent noise attenuation. Several classes of Radon transforms have been proposed. For instance, the linear Radon transform has been implemented to remove coherent noise [30] from seismic records. The parabolic Radon transform has been mainly utilized to separate multiples from primary reflections in Normal-Move-Out (NMO) corrected Common-Mid-Point (CMP) gathers [12]. Both linear and parabolic Radon transforms are usually implemented in the frequency-space (f − x) domain via the solution of an inverse problem that entails solving a linear system of equation for each frequency, independently [5].
Radon transforms that utilize hyperbolic integration paths have been also proposed to attenuate multiple reflections in CMP gathers and for seismic data reconstruction. Hyperbolic Radon transforms are time-variant transforms that cannot be implemented in the f − x domain. Due to the computational load of the hyperbolic Radon transform, seismic data are usually modified, either by an approximate NMO correction applied to the data [12] or by stretching the data along time axis [29], in order to used the parabolic Radon transform instead. Recently, a fast algorithm based on the butterfly algorithm has been proposed in [17] for generalized Radon transforms. A fast hyperbolic Radon transform has also been represented as convolutions in log-polar coordinates in [21]. [26] proposed to increase the resolution of the Radon transform by incorporating sparsity constraints. A similar idea was also introduced by [23] and [24] to improve the focusing power of the frequency domain Radon transforms. Recently, [10] presented deconvolutive Radon transform as a generalization of the conventional Radon transform and to increase the temporal resolution. Today, high resolution Radon transforms (often called sparse Radon transforms) are part of our arsenal of techniques for de-multiple and seismic data reconstruction [6,14,27,15]. In general, sparsity promoting iterative algorithms are utilized to estimate the Radon coefficients [27,19]. To a degree, the successful and efficient implementation of such algorithms relies on having access to two efficient canonical operators: the inverse Radon operator and its adjoint, which are needed to estimate the Radon domain coefficients that can model the data.
Traditional solvers for computing linear and parabolic Radon transforms rely on the Toeplitz structure of the time invariant Radon operator and solve the problem for each frequency independently [18,22]. In this way the solution of original large system is replaced by the solution of smaller systems [12], corresponding to each frequency. Therefore, it is implemented more efficiently than the time domain correspondence due to the scalability in parallel computing. However, in traditional f − x algorithms, the coefficients are regularized only in p direction.
The Fourier slice theorem [20] provides an efficient way for computing the linear Radon transforms without matrix inversion. However, to our knowledge, this theorem has not been used to compute more general time-invariant Radon transforms such as the parabolic Radon transform which is extensively utilized in the field of applied geophysics. Here, we generalize the Fourier slice theorem to make it applicable to any general time-invariant Radon transform. We derive an analytic expression between the Fourier coefficients of a function and the Fourier coefficients of its general Radon transform. In this way, very fast implementations of the inverse operator and its adjoint are possible via Fast Fourier Transform (FFT) algorithms. These two canonical operators are utilized for efficient implementation of the sparse Radon transform via split Bregman iterative method [11]. The optimal number of algorithm iterations is determined by the discrepancy principle, making the algorithm very efficient for processing large seismic data sets. Numerical results of noise attenuation and wavefield separation from synthetic and field seismic data show high performance of the proposed method.

1.1.
Outline. This work is structured as follows. In section 2 we briefly introduce the time-invariant Radon transform in continuous form. Then the Fourier slice theorem and its generalization are described in sections 2.1 and 2.2, respectively. An inverse transformation of the Radon transform by the generalized Fourier slice theorem is formulated in section 2.3. We extent our formulation for treating discrete data in section 3 and analyse the complexity of the algorithm in section 3.1. In section 4 we solve sparse Radon transform by using a Bregman iterative method and the generalized Fourier slice theorem. We give synthetic examples for linear and parabolic Radon transforms and the stability against noise in section 5.1 and compare computational cost and recovery performance with traditional direct methods. An example of spatial interpolation of irregular data is presented in section 5.2. Finally, the method is applied to field data for wavefield separation in section 5.3 and conclusions are given in section 6.
2. The continuous time-invariant Radon transform. Through this paper, we use the terminology adopted in seismic data processing applications. Let v(t, x) be a seismogram (function) of spatial variable (offset) x and time t. The variable x is used to indicate the distance between a source and receiver in a variety of acquisition scenarios. The time-invariant Radon transform of v can be defined by the following integral where τ is intercept time, p is a parameter, and function ϕ determines the integration path. Specific seismic events contained in the seismogram v can be modelled by a family of parametric curves τ + pϕ(x) for a predefined function ϕ. For example, linear coherent noise and surface waves can be modelled by the linear Radon transform: ϕ(x) = x. Similarly, reflections (after partial NMO correction) can be modelled by parabolic curves: ϕ(x) = x 2 , [12]. Therefore, some events contained in the seismogram could be focused by (1) in the transform domain. The focused events in the Radon domain can be easily isolated and filtered. Formula (1) is a generalization of the slant stacking technique used in geophysical community. Considering the wide application of Radon transform defined in (1) efficient and accurate implementation of it has always been a goal for geophysicists. Current algorithms solve time-invariant Radon transform (1) via an inverse strategy in the f − x domain as which is solved for each frequency f , separately. Therefore, in the case of discrete data, one needs to solve a system of complex linear equations in order to recover a frequency slice of the Radon plane [5,12,18,25,30,22]: where v(f ) and u(f ) denote frequency slices of data and Radon panels, and L(f ) is a frequency-dependent linear operator. Due to variation of the operator with frequency, processing all frequencies can be a time-consuming task which can limit the usage of the transform for large data sets. Here, we call this method direct Radon transform and use it as a benchmark for comparing the performance of the new method presented in this paper. The Fourier slice theorem is an efficient tool for fast computation of Radon transform when ϕ(x) = x. Here, we generalize this theorem for fast computation of a general form of Radon transform.
2.1. The Fourier-slice theorem. In the case of linear Radon transform, the following theorem provides a simple way for computing the f − p coefficients of the Radon transform. Proof. The proof is simple and can be found in many textbooks (e.g, [16]).
This theorem establishes the relation (4) k x = f p between wave number k x , frequency f , and slowness p and states that the Fourier transformation of the projection of v along a slope p is equal to the radial slice taken through the 2-D Fourier transform of v with the same slope (see Figure 1). A fast and invertible algorithm for the linear Radon transform based on this theorem and pseudo-polar Fourier transform is given in [1,2]. A fast algorithm has also been proposed for 3-D Radon transform via backprojection in [3]. However, this theorem can not be utilized for computing the Radon transform (1) for a general form of ϕ.
In the next section, we generalize the theorem for this task.
v(t, x) Figure 1. Graphical representation of the Fourier slice theorem.
The seismogram v includes a linear event. The projection of v along the slope of the event is the red line shown in u which is the IFT of the radial slice taken fromv using relation k x = f p.

Generalized Fourier-slice theorem.
Theorem 2.2 (generalized Fourier-slice theorem). Let v(t, x) be a function having 1-D Fourier transformv(f, x), and u(τ, p) be its general Radon transform defined in (1) with 2-D Fourier transformû(f, k p ), then there is the simple relation between k p , f , and x.
Proof. Applying the 2-D Fourier transform to both sides of (1) giveŝ which confirms the relation (5).
According to the new theorem the Fourier transformation of a trace at offset x is the same as the radical slice through the 2-D Fourier transformation of the Radon transform with the slope ϕ(x). Figure 2 schematically shows the relation between computation of the linear Radon transform via (1) and via the generalized Fourierslice theorem. The difference between the Fourier-slice theorem and its generalized form is that the former relates the 2-D Fourier transformation of the data to the 1-D Fourier transformation of its Radon plane while the latter relates the 1-D Fourier transformation of the data to 2-D Fourier transformation of its Radon plane.

Inverse transform. The inverse Fourier transform recovers
and using the following expression we arrive to the following results = û(f, p)e i2πf τ df dp = u(τ = t − pϕ(x), p)dp .
Inverse Problems and Imaging Volume 11, No. 3 (2017), 501-519 3. The discrete time-invariant Radon transform. The discrete time-invariant Radon transforms are of particular interest in practical seismic applications when dealing with discrete data in time and space. We assume that the sampling of the data domain and the Radon domain are uniform. We designate ∆x, ∆t, ∆p, and ∆τ the sampling interval along x, t, p, and τ axes, respectively. We also define n x , n t , n p , and n τ the corresponding discrete index parameters, then the correspondance between the discrete and continuous variables is given by Usually, the lower limits t min , x min and τ min are zero. However, parameter p min can get negative values corresponding to, e.g., negative slopes or curvatures. We use 1 and 2 to define uniformly sampled frequency and wavenumber of the Radon plane: f ( 1 ) = 1 ∆f and k p ( 2 ) = 2 ∆k p , for 1 = 0, · · · , N τ − 1 and 2 = 0, · · · , N p − 1. Therefore, according to theorem 2.2 Having computedû the Radon transform u is obtained by 2-D inverse Fourier transformation. Figure 3 shows the points ofû(f, k p ) over imposed to the Cartesian grids ofv(f, x) for the linear and parabolic Radon transforms. Clearly, the computation ofû from the samples ofv requires an interpolation stage. The latter is implemented in the frequency domain and discussed with more detail as we progress with our exposition. The inverse Radon transform is performed by computingv(f, x) from u(f, k p ) via the expression followed by inverse Fourier transform along the frequency axis. Let f max be the maximum frequency in the trace at largest offset, x max , then according to the relation (5) the maximum value of k p is defined as k pmax = f max ϕ(x max ), and thus which is in accordance with the result of [25]. Furthermore, since k pmax = (N p − 1)∆k p we obtain (19) ∆k On the other hand, according to the Nyquist sampling theory the maximum value of p recoverable from the sampled spectrum is defined as and aliasing occurs for the events with p values above this. For the computed spectrum of length N p we need to perform an inverse DFT of length, at least, N p and naturally the resulting Radon coefficients will be N p -periodic. The left and right halves of the output sequence are thus swapped to get the Radon coefficients corresponding to the range [−p max , p max ].  N p logN τ N p ), in order to get the estimated Radon coefficients. The same mechanism with reversed order can also be used for the inverse Radon transform to reconstruct the data from the Radon coefficients. Therefore, the implementation of both the forward and inverse Radon transformations has a complexity given by O (N logN ) where N is the maximum number of samples in the data space or model space. Generally, the number of samples in the data space is larger than or equal to the number of samples in the Radon space. Thus, N can be considered as the number of data samples. 4. Sparse Radon transforms. Time-invariant Radon transforms have been extensively used for noise attenuation and wavefield separation. The method is based on the focusing power of the Radon transform [27]. However, a practical implementation of the transform requires to consider the finiteness of the data and discretization effects that often deteriorate the resolution of the resulting Radon coefficients. Therefore, regularization strategies are required in order to estimate highly resolved (focused) distribution of the Radon coefficients. Using the notation from linear algebra and based on the steps of the algorithm discussed above, the inverse Radon transform via the generalized Fourier slice theorem can be written as (22) v = F −1 1 SF 2 u where F 1 and F 2 denote 1-D FFT and 2-D FFT operators of appropriate size. The operator S re-samples the Fourier coefficients of the Radon transform from the Cartesian grid to the grid points defined by (17). Figure 3 shows the relation between the gridpoints ofû and those ofv for forward/inverse linear and parabolic Radon transforms. The figure is shown for only the positive frequencies and wavenumbers. In the case of the linear transformation, for each frequency f , the distance between successive points along the wavenumber axis remains constant. These points form a so called pseudopolar grid and [1,2] used the chirp z-transform to calculate the coefficients value at these points. In the case of parabolic Radon transform, however, the distance between successive points are non-uniform and demands another interpolation scheme. As seen from Figure 3, in the inverse transformation, the first column and the diagonal coefficients ofû correspond to the first and last traces, respectively, and thus there is no need for interpolation. One needs to interpolate the coefficients of the intermediate traces. The algorithm in this paper employs linear interpolation for its efficiency and acceptable accuracy. The more smoother the data under analysis is the better will be the performance of linear interpolation. Generally, seismic data are band-limited both in time and in spatial directions. Therefore, the method works satisfactorily in most situations of seismic studies. However, one may use more accurate interpolation methods for complicated data.
Using the formulation given by (22), the Radon coefficients can be regularized in order to increase its resolution by using an appropriate sparsity promoting regularization. Introducing a new variableb = F 1 v (which is indeed the f − x domain representation of the data) reduces problem (22) to a simple interpolation problem Sû =b, whereû = F 2 u. Therefore, the problem of finding a sparse solution can be defined as u = F −1 2û whereû is a solution to the basis pursuit program [8] arg min u ||F −1 2û || 1 subject to Sû =b. (24) where || · || 1 denotes the sum of absolute values. By this program we look for the sparsest Radon panel, in the one-norm sense, whose 2-D FFT at the specified grid points (see Figures 3b and d) matches the f −x coefficients of the data,b. However, the presence of uncertainties in the data or inaccuracy in the interpolation may causê b out of range and hence (24) may not have a solution. A better choice is arg min where ε is the error bound. For linear interpolation, the matrix of S is sparse (there will be two non-zero entries at each row) and when combined with efficient applications of F 2 /inverse F 2 via FFT's make efficient solution of problem (25) possible via split Bregman iterations [11]. The details of the algorithm is fully analysed in [11]. We just provide the algorithm (Algorithm 1) for the solution of problem (25). In the algorithm shrink is the soft shrinkage function and d is the diagonal elements of αS T S + βI. The advantage of the split Bregman iterations employed for solving our problem is that αS T S +βI can be well approximated by its diagonal, where I denotes identity matrix. Therefore, the run time for application of each Bregman iteration is mostly due to 2-D FFT/IFFT's and hence the total cost at each iteration is O (N logN ) where N is the number of coefficients. Aside from fast convergence of the Bregman iteration a main advantage of it for solving our problem is that here the number of Bregman iterations has the role of regularization parameter. We allow the iterations to proceed until the misfit between predicted and observed data reach a desired predefined tolerance.   5.1. Synthetic data. Two data sets have been simulated in order to examine the functionality of the proposed Radon transforms. The first (Figure 4a) includes linear events and the second (Figure 4b) consists of parabolic events. Both data sets are of dimension 1024 by 1024. The absolute value of the f − x representation of the data are shown below each data set. We analyzed the data for f max = 50 Hz and the resulting paths of Radon coefficients are shown in green in Figure 4. Both non-sparse and sparse Radon transforms have been applied to the simulated data using the parameters defined above and the results are shown in Figures 5 and  6 for the linear and parabolic transforms, respectively. The predicted data from the estimated Radon transform is also shown below each figure. It can be seen from our results that for both non-sparse and sparse methods the estimated Radon  coefficients can predict the data with a high degree of accuracy. The non-sparse coefficients, as expected, clearly show energy leakage from the main peaks along two crossing paths. On the other hand, the sparse coefficients are focused and free of smearing. From a computational time viewpoint, however, the non-sparse coefficients have been obtained in 28e-3 seconds while the sparse coefficients have been calculated by 20 Bregman iterations with a total time of 1.2 sec. We continued by applying the parabolic Radon transform to data of different sizes. Both forward (sparse and non-sparse) and inverse transforms were examined. In order to perform the new method, we employed linear interpolation to calculate the coefficients at desired points as in Figure 3. To show the effect of problem size, all data sets examined were square, N t = N x , with equal number of samples in the data and model space and the analysis was done for all frequencies up to the Nyquist. We compared the new method with the traditional f − x algorithm to evaluate the least-squares forward transform and employed three different fast linear solvers to solve the corresponding complex Toeplitz systems in 3 for each frequency: the MATLAB backslash operator, the Levinson recursion, and the preconditioned conjugate gradient method with the Chan preconditioner (PCG-Chan) [7]. A constant value has been added to the main diagonal of the Toeplitz matrices to stabilize the inversion in the forward transform. Table 1 shows computational times, speed-up of the proposed method over these three methods, and the meansquared-error of the results of the proposed method with respect to the results of  the backslash method. All results shown were averaged over 20 runs. One can see that the new method accurately performed significantly faster than the traditional fast linear solvers. The data has also been transformed by sparse parabolic Radon transform using Algorithm 1. Table 2 shows the number of Bregman iterations, average time per iteration, and total computational time (in sec) required for convergence to a given misfit of 10 −3 . Here again all results shown were averaged over 20 runs. It is seen that the number of Bregman iteration does not increase with the size of data. Furthermore, a comparison of the run times in Table 2 with those reported in Table  1 shows that the proposed algorithm for sparse transform is even faster than the traditional non-sparse solvers.
The new method has also been tested for inverse parabolic Radon transform and compared with the traditional direct calculations defined in 3. The results are shown in Table 3 for different sizes of data. The inverse transform has been calculated significantly faster by the generalized Fourier slice theorem.
We continued by performing the new method in the presence of noise. The gathers ( Figure 4) were contaminated by additive random noise of signal-to-noiseratio (SNR) -5 dB (Figures 7a and b). Here again we analysed the data for f max =

5.2.
Interpolation of irregular data. Seismic data sets are generally irregularly sampled in the spatial direction due to variable operating conditions during seismic data acquisition. This irregular sampling can limit the effectiveness of data processing and imaging algorithms. In general, acquired traces need to be interpolated to a regular grid prior to start classical processes for signal enhancement and data inversion [9]. The Radon transform is a useful and efficient tool for seismic data regularization. The sampling operator S in Algorithm 1 explicitly depends upon the spacial coordinates of the traces x, therefore, when the Radon coefficients have been estimated, the interpolated data can easily be constructed by the inverse Radon transformation to any desired grid points. We use 24 traces of the data shown in Figure 4 collected at random locations (see Figures 8a and b) and used Algorithm 1 to perform sparse Radon transforms. Then, the inverse transformation was used to synthesize data at all grid points. The original data before and after the process of reconstruction are shown in Figure 8. It can be seen that the original data sets have been accurately constructed.
5.3. Field data. As a final example, we tested our method on a real CMP gather and compared it with the direct method. The parabolic Radon transform is known as a successful tool for the attenuation of multiple reflections in CMP gathers [12]. Figure 9 shows a field CMP gather after NMO correction applied using the velocity of primary reflections. As can be seen from this figure primary reflections are flatted but multiple reflections show parabolic moveout. Figure 10a shows the sparse parabolic Radon panel generated by the proposed method (300 iterations of Algorithm 1) and Figure 10b shows the corresponding predicted data. As can be seen, the energy of primary reflections are well clustered around p value zero and separated form the multiples energy, enclosed by the rectangle in Figure 10a. We just transformed back the multiples energy ( Figure 10c) and subtracted it from the original data to get the primary only wavefield as shown in Figure 10d. The same procedure was applied to the data using traditional direct method and 300 iterations of fast iterative shrinkage-thresholding algorithm (FISTA) [4] and the results are shown in Figure 11. As can be seen from Figures 10 and 11, the obtained results by the two methods are very similar. However, the average time per iteration for the proposed method was 0.01 sec while it was 0.6 sec for direct method.