A fast butterfly algorithm for generalized Radon transforms

Introduction

In seismic data processing, the Radon transform (RT) (Radon, 1917), or slant stack, is a set of line integrals that map mixed and overlapping events in seismic gathers to a new transformed domain where they can be separated (Gardner and Lu, 1991). The integrals can also be taken along curves: parabolas (parabolic RT), or hyperbolas (hyperbolic RT or velocity stack) are most commonly used. A major difference between these transforms is that the former two are time-invariant (i.e., involve a convolution in time) whereas the latter is time-variant. When the curves are time-invariant, the transform can be performed efficiently in the frequency domain using the convolution theorem. However, this approach does not work for time-variant transforms. As a result, the hyperbolic Radon transform is usually thought of as requiring a computation in the time domain, which is computationally expensive due to the large size of seismic data. Nevertheless, the hyperbolic transform is often preferred as it better matches the true seismic events in common midpoint (CMP) gathers (Thorson and Claerbout, 1985).

In this work, we construct a fast butterfly algorithm to effectively evaluate time-variant transforms such as the hyperbolic Radon transform. As opposed to the conventional, relatively costly velocity scan (i.e., direct integration + interpolation in the time domain), our method provides an accurate approximation in only $O(N^2\log N)$ (all the $\log$ in this paper refer to logarithm to base 2) operations for 2D data. Here $N$ depends on the maximum frequency and offset in the dataset and the range of parameters (intercept time and slowness) in the model space, and can often be chosen small compared to the grid size. The adjoint of the transform can be evaluated similarly without extra difficulty. Note that the algorithm introduced in this paper only deals with the fast implementation of a single integral operator (forward Radon transform or its adjoint), not an iteration process for its inversion which is the main objective of many previous works on fast Radon transforms (Liu and Sacchi, 2002; Wang and Ng, 2009; Trad et al., 2002; Sacchi, 1996).

Radon transforms have been widely used to separate and attenuate multiple reflections (Moore and Kostov, 2002; Hargreaves et al., 2003; Herrmann et al., 2000; Trad, 2003; Foster and Mosher, 1992; Hampson, 1986; Yilmaz, 1989). As having fast implementations of both forward and adjoint transforms is an essential component of least-squares minimization, our hope is that the current fast algorithm will help to make the hyperbolic Radon transform an accessible tool for improving the inversion process.

The term ``generalized Radon transform" connotes a broader context where integrals are taken along arbitrary parametrized sets of smooth curves. The term was introduced by Beylkin (1984,1985), who showed that an asymptotically-correct inverse follows from an amplitude correction to the adjoint. Kirchhoff migration and its (regularized) inverse can be expressed as generalized Radon transforms. The algorithm presented in this paper can in principle be applied in the context of Kirchhoff migration, although we do not attempt to do so here.

The rest of the paper is organized as follows. We first introduce the low-rank approximations and the butterfly structure of the hyperbolic Radon operator; then use these building elements to construct our fast algorithm. A brief description of the algorithm is given in the main text, and a complete derivation can be found in the appendix. We present numerical examples using both synthetic and field data to illustrate the accuracy and efficiency of the proposed algorithm.

A fast butterfly algorithm for generalized Radon transforms