Fast hyperbolic Radon transform represented as convolutions in log-polar coordinates

doi:10.1016/j.cageo.2017.04.013

Computers & Geosciences

Volume 105, August 2017, Pages 21-33

https://doi.org/10.1016/j.cageo.2017.04.013 Get rights and content

Highlights

•
Fast GPU implementation of the log-polar-based method for computing hyperbolic Radon transforms.
•
Multiple attenuation and interpolation using fast hyperbolic Radon transforms.
•
Accuracy and performance tests for the log-polar-based method of computing hyperbolic Radon transforms.

Abstract

The hyperbolic Radon transform is a commonly used tool in seismic processing, for instance in seismic velocity analysis, data interpolation and for multiple removal. A direct implementation by summation of traces with different moveouts is computationally expensive for large data sets. In this paper we present a new method for fast computation of the hyperbolic Radon transforms. It is based on using a log-polar sampling with which the main computational parts reduce to computing convolutions. This allows for fast implementations by means of FFT. In addition to the FFT operations, interpolation procedures are required for switching between coordinates in the time-offset; Radon; and log-polar domains. Graphical Processor Units (GPUs) are suitable to use as a computational platform for this purpose, due to the hardware supported interpolation routines as well as optimized routines for FFT. Performance tests show large speed-ups of the proposed algorithm. Hence, it is suitable to use in iterative methods, and we provide examples for data interpolation and multiple removal using this approach.

Introduction

In the processing of Common-Midpoint gathers (CMPs), the hyperbolic Radon transform has proven to be a valuable tool for instance in velocity analysis (Clayton and McMechan, 1981, Greenhalgh et al., 1990); aliasing and noise removal (Turner, 1990); trace interpolation (Averbuch et al., 2001, Yu et al., 2007); and attenuation of multiple reflections (Hampson, 1986). The hyperbolic Radon transform is defined as $R_{h} f (τ, q) = \int_{- \infty}^{\infty} f (\sqrt{τ^{2} + q^{2} x^{2}}, x) dx,$ where the function $f (t, x)$ usually corresponds to a CMP gather. Here, the parameter q characterizes an effective velocity value; and τ represents the intercept time at zero offset.

Several versions of Radon transforms are used in seismic processing, e.g., straight-line, parabolic, and hyperbolic Radon transforms. In many applications there is a need for a sparse representation of seismic data using hyperbolic wave events. One way to get sparse representations is by using iterative thresholding algorithms with sparsity constraints (Daubechies et al., 2004, Sacchi and Ulrych, 1995). Popular applications using such representations are seismic data interpolation and wavefield separation (Jiang et al., 2016, Trad, 2003). Since iterative schemes for computing such representations require the application of the forward and adjoint operators several times, it becomes important to use fast algorithms to limit the total computational cost.

Note that the direct summation over hyperbolas in (1) has a computational complexity of $O (N^{3})$ , given that the numbers of samples for the variables $t, x, τ, q$ are $O (N)$ . There are many effective ( $O (N^{2} \log N)$ ) methods for rapid evaluation of the traditional Radon transforms, or the parabolic Radon transform, see Beylkin, 1984, Fessler and Sutton, 2003, Schonewille and Duijndam, 2001. The hyperbolic Radon transform is, however, more challenging. Nonetheless, a fast ( $O (N^{2} \log N)$ ) method for hyperbolic Radon transforms was recently presented in Hu et al. (2013). The method is based on using the fast butterfly algorithms described in O'N, 2007, Candès et al., 2009, and versions addressing computational efficiency are presented in Poulson et al., 2014, Li et al., 2015b, Li et al., 2015a, Li and Yang, 2016.

A fast method for the standard Radon transform was proposed in Andersson (2005) by expressing the Radon transform and its adjoint in terms of convolutions in log-polar coordinates. To use this approach, it is necessary to resample data in log-polar coordinates, and this requires some interpolation method. An important property of such a scheme is that since the interpolation procedures are performed in the Radon and image domains, the errors will be kept local. On the contrary, the interpolation errors in the frequency domain will create global non-systematic errors when switching back to the time-space domain. Since these non-localized errors contribute to the general structure of the function, they have to be mainly controlled.

Computationally efficient algorithms for GPUs were presented for the log-polar-based approach in Andersson et al. (2016). In this paper we propose to use the same approach and construct algorithms with complexity $O (N^{2} \log N)$ for evaluation of the hyperbolic Radon transform. We present computational performance tests confirming the expected accuracy and the computational complexity, as well as predicted computational speed-ups for parallel implementations. Finally, we present several synthetic and real data tests using the hyperbolic Radon transform for data interpolation and multiple attenuation.

Section snippets

Method

To begin with, we note that functions $f (t, x)$ describing CMP gathers are symmetric with respect to x=0. Hence, by introducing $\tilde{f} (s, y) = \frac{f (\sqrt{s}, \sqrt{y})}{2 \sqrt{y}},$ it follows that $R_{h} f (τ, q) = 2 \int_{0}^{\infty} f (\sqrt{τ^{2} + q^{2} x^{2}}, x) dx = 2 \int_{0}^{\infty} \tilde{f} (τ^{2} + q^{2} y, y) dy .$ The resulting expression in (3) has a form of the Radon transform over straight lines, and a fast algorithm for the evaluation of this was presented in Andersson et al. (2016), referred to as the log-polar Radon transform which is based on rewriting the key operations as convolutions in a

Reconstruction techniques

The adjoint operator for the hyperbolic Radon transform $R_{h}^{*}$ is defined by using the inner product equality $〈 R_{h} f, g 〉 = 〈 f, R_{h}^{*} g 〉,$ for arbitrary f and g. The operator is easy to construct by using the approach with switching to log-polar coordinates, essentially by reversing the order of the operations. With the adjoint operations at hand, one can consider iterative methods for representing f by sparse sums of hyperbolic wave events, and related interpolation and reconstruction techniques. A popular

Discretization

In this section we derive guidelines for how to choose discretization parameters. For simplicity, we assume to work with regular sampling in the time-offset, and in the Radon domain; but the log-polar-based method can be easily generalized for unequally-spaced grids in these two domains.

In order to apply FFTs, samples in log-polar coordinates $(θ, ρ) \in [- \frac{β}{2}, \frac{β}{2}] \times [\log a_{r}, 0]$ must be chosen on an equally spaced grid. By using coordinate conversions for the log-polar setup given by $(\begin{matrix} φ (t, x) \\ η (t, x) \end{matrix}) = P_{1} T (\begin{matrix} t^{2} \\ x^{2} \end{matrix})$

Accuracy tests

The proposed method uses cubic B-splines for data interpolation and therefore does not have provable error estimates. Here we aim at seismic applications in which case the data are band-limited (in particular, band-pass filtering is applied at the preprocessing step) and a limited accuracy level of several digits is sufficient in most applications. Thus in this section we study accuracy levels when processing synthetic and real CMP gathers which have different frequency content.

To begin with,

Performance tests

In this section we will check the performance of the log-polar-based method and compare it to the performance of the fast butterfly algorithm. The tables in the previous section show that the log-polar-based method using the cubic interpolation and one splitting in variables (according to (18)) behave similarly in terms of accuracy as the fast butterfly algorithm with q_i=9 Chebyshev points and $M / 16$ as the maximum value of the phase function. Thus, we compare the computational times for these

Applications

In this section we mention some applications of the fast hyperbolic Radon transform. These are fairly standard, but the examples could be of practical interest due to the substantial computational speedup of the proposed implementation of the hyperbolic Radon transform.

Non-hyperbolic moveouts

In this work we have dealt with the hyperbolic moveout $t (x) = \sqrt{τ^{2} + q^{2} x^{2}}$ . The proposed method can be trivially adapted for the parabolic Radon transform where the moveout is given by $t (x) = τ + {qx}^{2}$ (can be used for small offsets or after the normal moveout correction). Several other moveout types are used in seismic processing as listed in Fomel and Stovas (2010) . Let us mention a few examples of these non-hyperbolic moveouts: $t {(x)}^{2} = q_{1} + q_{2} x^{2} + q_{3} x^{4}$ and $t {(x)}^{2} = q_{1} + q_{2} x^{2} + q_{3} x^{4} / (1 + q_{4} x^{2})$ (associated with seismic

Conclusions

A fast log-polar-based method for the evaluation of the hyperbolic Radon transform has been presented. According to the tests performed, the method demonstrates reasonable accuracy and favorable computational costs compared to other methods. The accuracy of the method can be increased when considering higher order interpolation kernels for coordinate conversions between time-offset, Radon, and log-polar domains. Numerical tests show that the GPU implementation is more than 10,000 faster for

Acknowledgements

This work was supported by the Crafoord Foundation (2014-0633) and the Swedish Research Council (2011-5589 and 2015-03780)

References (35)

A. Averbuch et al.
Fast and accurate polar Fourier transform
Appl. Comput. Harmon. Anal.
(2006)
G. Beylkin
On the fast Fourier transform of functions with singularities
Appl. Comput. Harmon. Anal.
(1995)
S. Greenhalgh et al.
Seismic wavefield separation by multicomponent tau-p polarisation filtering
Tectonophysics
(1990)
H. Lu et al.
Three-dimensional architecture of shelf-building sediment drifts in the offshore canterbury basin, New Zealand
Mar. Geol.
(2003)
F. Andersson
Fast inversion of the Radon transform using log-polar coordinates and partial back-projections
SIAM J. Appl. Math.
(2005)
F. Andersson et al.
Fast algorithms and efficient GPU implementations for the Radon transform and the back-projection operator represented as convolution operators
SIAM J. Imaging Sci.
(2016)
A. Averbuch et al.
Fast slant stack: a notion of radon transform for data in a cartesian grid which is rapidly computible, algebraically exact, geometrically faithful and invertible
SIAM J. Sci. Comput.
(2001)
S. Basu et al.
$O (N^{2} \log N)$ filtered backprojection reconstruction algorithm for tomography
IEEE Trans. Image Process.
(2000)
G. Beylkin
The inversion problem and applications of the generalized Radon transform
Commun. Pure Appl. Math.
(1984)
E. Candès et al.
A fast butterfly algorithm for the computation of Fourier integral operators
Multiscale Model. Simul.
(2009)

R.W. Clayton et al.

Inversion of refraction data by wave field continuation

Geophysics

(1981)

I. Daubechies et al.

An iterative thresholding algorithm for linear inverse problems with a sparsity constraint

Commun. Pure Appl. Math.

(2004)

C. De Boor

A Practical Guide to Splines

(1978)

A. Dutt et al.

Fast Fourier transforms for nonequispaced data

SIAM J. Sci. Comput.

(1993)

J.A. Fessler et al.

Nonuniform fast fourier transforms using min-max interpolation

IEEE Trans. Signal Process.

(2003)

S. Fomel et al.

Generalized nonhyperbolic moveout approximation

Geophysics

(2010)

S. Fomel et al.

Madagascar: open-source software project for multidimensional data analysis and reproducible computational experiments

J. Open Res. Softw.

(2013)

Cited by (19)

A local radon transform for seismic random noise attenuation
2021, Journal of Applied Geophysics
Citation Excerpt :
What is more, a new fast Radon transform method which is based on compressive sensing, a nonlinear seismic sampling technique, has been used in producing high-resolution results consuming less computational time and the number of measurements (Latif and Mousa 2016). Furthermore, it is proved that the Graphical Processor Units (GPUs) has its spacial function in dealing with the statistics with big data and complex computational processing, which plays an important role in the velocity analysis and also multiple attenuation (Nikitin et al. 2017). Sacchi and Ulrych (1995) realized sparse RT in the frequency domain by using the probability-density function of Cauchy form and Bayesian rule, which led to an operator that does not exhibit a Toeplitz structure.
The Radon transform can sparsify the seismic data assuming that the events follow a certain type of trajectories. Commonly known Radon transforms include linear, parabolic, and hyperbolic transforms, which are usually exploited to process pre-stack data with the parabolic and hyperbolic kernels, or structurally simple post-stack seismic data with distinct linear events. We develop a new Radon transform by introducing a local sparsification strategy into the traditional method. We first select a set of spatial positions, where the local Radon transform is centered in, then we construct the Radon transform with a fixed width; the linear-kernel is shared by all local Radon transforms. The local Radon coefficients are obtained directly from solving an augmented Radon transform model. To enforce the sparsity of coefficients, we apply an L1-norm regularization to the inverse model and solve it via a preconditioned least-squares method. We use a seismic shot gather with curved events and a complicated post-stack field dataset to demonstrate the superiority of the proposed method in dealing with complicated seismic data compared with the conventional Radon transform.
Dynamic in-situ imaging of methane hydrate formation and self-preservation in porous media
2020, Marine and Petroleum Geology
Citation Excerpt :
The constructed pipeline for studying gas-hydrate formation processes includes plugins for flat field correction, ring removal, and phase-retrieval filter. We also prepared our own plugin for tomography data reconstruction by using the fast log-polar-based imaging algorithm optimized for GPUs (Andersson et al., 2016; Nikitin et al., 2017). Following Lei et al. (2018) we conducted experiments in order to choose optimal water salinity for better contrasts between different phases: sand grains, water, methane hydrate, gas.
We present the results of dynamic in-situ 3D X-ray imaging of methane hydrates microstructure during methane hydrate formation and dissociation in sand samples. Short scanning times and high resolution provided by synchrotron X-rays allowed for better understanding of water movement and different types of gas-hydrate formation. Complementing previous observations, we conclude that the process of gas-hydrate formation is accompanied by the water movements caused by cryogenic water suction that happens in sequences of short fast movements with longer equilibrium states in between (when the water is immobile). Based on the 3D microstructure we identified two distinct types of gas-hydrate formation: (i) into the gas pockets and (ii) inside water volumes. For both mechanisms we do not see problems in gas or water supply to support the gas-hydrate formation. The rate of dissociation in the self-preservation mode (pressure drop at negative temperatures) appears to be different for these two types of gas hydrates. This means that the history of the gas-hydrate formation may influence its behaviour at the dissociation stage (e.g. gas-hydrate production).
Fast hyperbolic Radon transform using chirp-z transform
2019, Digital Signal Processing: A Review Journal
The hyperbolic Radon transform plays an important role in seismic data processing for its ability to focus seismic events in the transform domain. Traditional algorithms based on direct implementations, however, are inefficient with limited applications for processing large data sets. A new algorithm is presented for fast computation of the hyperbolic Radon transform and its sparse calculations. It uses interpolation procedures to stretch the data along time axis and efficiently computes the summation paths in the new coordinates via the chirp-z transform which is carried out by fast Fourier transform (FFT). The proposed fast algorithm is then used within the deconvolutive form of the Radon transform and iterative sparse algorithms for effective decomposition of CMP gathers with an improved temporal resolution, compared to the traditional Radon transforms. The effectiveness of the new algorithm are confirmed on sparse velocity-stack inversion, primary and multiple separation, high-quality stacking, and automatic velocity model building. The tests show that sparse velocity-stack inversion using the new algorithm is even more efficient than the traditional velocity scan, both in resolution and speed. Furthermore, numerical tests show the superiority of the proposed algorithm over the state-of-the-art fast algorithms, based on butterfly scheme and log-polar convolutions, demanding less computational complexity.
A geometric partitioning method for distributed tomographic reconstruction
2019, Parallel Computing
Citation Excerpt :
We distinguish between two approaches that are being taken in algorithm research for fast tomography. First, alternative reconstruction algorithms are being developed that approximate advanced but slow iterative methods, by faster and lighter methods [15–19]. Second, techniques are being developed that take advantage of advances in computer hardware.
Tomography is a powerful technique for 3D imaging of the interior of an object. With the growing sizes of typical tomographic data sets, the computational requirements for algorithms in tomography are rapidly increasing. Parallel and distributed-memory methods for tomographic reconstruction are therefore becoming increasingly common. An underexposed aspect is the effect of the data distribution on the performance of distributed-memory reconstruction algorithms. In this work, we introduce a geometric partitioning method, which takes into account the acquisition geometry and aims to minimize the necessary communication between nodes for distributed-memory forward projection and back projection operations. These operations are crucial subroutines for an important class of reconstruction methods. We show that the choice of data distribution has a significant impact on the runtime of these methods. With our novel partitioning method we reduce the communication volume drastically compared to straightforward distributions, by up to 90% for a number of cases, and furthermore we guarantee a specified load balance.
The importance of including density in multiparameter asymptotic linearized direct waveform inversion: A case study from the Eastern Nankai Trough
2022, Geophysical Journal International
A Fast Sparse Hyperbolic Radon Transform Based on Convolutional Neural Network and Its Demultiple Application
2022, IEEE Geoscience and Remote Sensing Letters

View all citing articles on Scopus

View full text

Research paperFast hyperbolic Radon transform represented as convolutions in log-polar coordinates

Highlights

Abstract

Introduction

Section snippets

Method

Reconstruction techniques

Discretization

Accuracy tests

Performance tests

Applications

Non-hyperbolic moveouts

Conclusions

Acknowledgements

Appl. Comput. Harmon. Anal.

Appl. Comput. Harmon. Anal.

Tectonophysics

Mar. Geol.

Fast inversion of the Radon transform using log-polar coordinates and partial back-projections

SIAM J. Appl. Math.

Fast algorithms and efficient GPU implementations for the Radon transform and the back-projection operator represented as convolution operators

SIAM J. Imaging Sci.

Fast slant stack: a notion of radon transform for data in a cartesian grid which is rapidly computible, algebraically exact, geometrically faithful and invertible

SIAM J. Sci. Comput.

O(N2logN) filtered backprojection reconstruction algorithm for tomography

IEEE Trans. Image Process.

The inversion problem and applications of the generalized Radon transform

Commun. Pure Appl. Math.

A fast butterfly algorithm for the computation of Fourier integral operators

Multiscale Model. Simul.

Inversion of refraction data by wave field continuation

Geophysics

An iterative thresholding algorithm for linear inverse problems with a sparsity constraint

Commun. Pure Appl. Math.

A Practical Guide to Splines

Fast Fourier transforms for nonequispaced data

SIAM J. Sci. Comput.

Nonuniform fast fourier transforms using min-max interpolation

IEEE Trans. Signal Process.

Generalized nonhyperbolic moveout approximation

Geophysics

Madagascar: open-source software project for multidimensional data analysis and reproducible computational experiments

J. Open Res. Softw.

Research paper
Fast hyperbolic Radon transform represented as convolutions in log-polar coordinates

$O (N^{2} \log N)$ filtered backprojection reconstruction algorithm for tomography