next up previous [pdf]

Next: Low-rank approximations Up: A fast butterfly algorithm Previous: Introduction

Algorithm

Assume $ d(t,h)$ is a function in the data space. The hyperbolic Radon transform $ R$ maps $ d$ to a function $ (Rd)(\tau,p)$ in the model space (Thorson and Claerbout, 1985) through

$\displaystyle (Rd)(\tau,p)=\int d(\sqrt{\tau^2+p^2h^2},h)\,dh,$ (1)

where $ t$ is time, $ h$ is offset, $ \tau$ is intercept time, and $ p$ is slowness. Fixing $ (\tau,p)$ , hyperbola $ t=\sqrt{\tau^2+p^2h^2}$ describes the traveltime for the event; hence integration along these curves can be used to identify different reflections.

Instead of approximating the integral in equation 1 directly, we reformulate it as a double integral,

$\displaystyle (Rd)(\tau,p)=\iint\hat{d}(f,h) e^{ 2\pi i f\sqrt{\tau^2+p^2h^2} } \,df\,dh.$ (2)

Here $ f$ is the frequency, $ \hat{d}(f,h)$ is the Fourier transform of $ d(t,h)$ in $ t$ . A simple discretization of equation 2 yields

$\displaystyle (Rd)(\tau,p)=\sum_{f,h} e^{ 2\pi i f\sqrt{\tau^2+p^2h^2} }\hat{d}(f,h)$ (3)

(the area element is omitted; the same symbols $ f$ , $ h$ , $ \tau$ , and $ p$ are used for both continuous and discrete variables). The reason that hyperbolic RT is harder to compute than linear RT ($ t=\tau+ph$ ) or parabolic RT ( $ t=\tau+ph^2$ ) should be clear from equation 3: product $ f\tau$ in the phase cannot be decoupled from other terms.

To construct the fast algorithm, we first perform a linear transformation to map all discrete points in $ (f,h)$ and $ (\tau,p)$ domains to points in the unit square $ [0,1]\times[0,1]$ ( $ [a,b]\times[c,d]$ represents a 2D rectangular domain in the $ xy$ -plane with $ x\in[a,b]$ and $ y\in[c,d]$ ), i.e., a point $ (f,h)\in[f_$min$ ,f_$max$ ]\times[h_$min$ ,h_$max$ ]$ is mapped to $ \mathbf{k}=(k_1,k_2)\in[0,1]\times[0,1]=K$ via

$\displaystyle f=(f_$max$\displaystyle -f_$min$\displaystyle )k_1+f_$min$\displaystyle , \quad h=(h_$max$\displaystyle -h_$min$\displaystyle )k_2+h_$min$\displaystyle ;$ (4)

a point $ (\tau,p)\in[\tau_$min$ ,\tau_$max$ ]\times[p_$min$ ,p_$max$ ]$ is mapped to $ \mathbf{x}=(x_1,x_2)\in[0,1]\times[0,1]=X$ via

$\displaystyle \tau=(\tau_$max$\displaystyle -\tau_$min$\displaystyle )x_1+\tau_$min$\displaystyle , \quad p=(p_$max$\displaystyle -p_$min$\displaystyle )x_2+p_$min$\displaystyle .%\footnote{The notations $\mathbf{x}$ and $\mathbf{k}$ are adopted to be consistent with the Cand\\lq {e}s et al. paper referenced below.}
$ (5)

If we define input $ g(\mathbf{k})=\hat{d}(f(k_1),h(k_2))$ , output $ u(\mathbf{x})=(Rd)(\tau(x_1),p(x_2))$ , and phase function $ \Phi(\mathbf{x},\mathbf{k})=f(k_1)\sqrt{\tau(x_1)^2+p(x_2)^2h(k_2)^2}$ , then equation 3 can be written as

$\displaystyle u(\mathbf{x})=\sum_{\mathbf{k}\in K} e^{ 2\pi i \Phi(\mathbf{x},\mathbf{k})}g(\mathbf{k}), \quad \mathbf{x}\in X$ (6)

(throughout the paper, $ K$ and $ X$ will either be used for sets of discrete points or square domains containing them; the meaning should be clear from the context). This form is the discrete version of an oscillatory integral of the type

$\displaystyle u(\mathbf{x})=\int_K e^{ 2\pi i \Phi(\mathbf{x},\mathbf{k})}g(\mathbf{k}) \,d\mathbf{k},\quad \mathbf{x}\in X,$ (7)

whose fast evaluation has been considered in Candès et al. (2009). Our method for computing the summation in equation 6 follows the FIO butterfly algorithm introduced there.



Subsections
next up previous [pdf]

Next: Low-rank approximations Up: A fast butterfly algorithm Previous: Introduction

2013-07-26