The discrete Radon transform
In practice, a discrete form of the rho filter ρ(τ) can be convolved with the discrete form of u(v, τ) prior to summing over a finite range of velocities to reconstruct the original data d(h, t). Note that the rho filter is a 1-D filter that operates on each trace of u(v, τ), individually. It can account for the suppression of high frequencies during summation along discretely sampled traveltime curves, t = τ+ϕ(v, h). However, one can intuitively contend that the rho filtering alone cannot remove the amplitude smearing along the velocity axis in Figure 6.4-2d. Based on a rigorous mathematical analysis, Beylkin  showed that a discrete form of Radon’s inversion formula (equation 11b) can be expressed as a least-squares solution to a set of linear system of equations. We shall review the discrete radon transform for the special case of velocity-stack transformation.
To reduce the amplitude smearing on conventional velocity-stack gathers Thorson and Claerbout  proposed a least-squares formulation. There is one issue, however, we must take into consideration before we seek a practical scheme to compute the inverse Radon transform defined by equation (11b). To comply with the linear form the Radon transform defined by equation (11a), apply stretching in the time direction by setting t′ = t2 and τ′ = τ2 . Equation (9b) then takes the form
As a result of the t2-stretching, hyperbolic events in the offset domain are transformed to parabolic events described by the traveltime equation (12). Note that, unlike the hyperbolic events, the moveout 4h2/v2 associated with the parabolic events is time-independent.
The objective is to estimate the transform u(v, τ′) such that the difference between the actual CMP gather d(h, t′) and the modeled CMP gather d′(h, t′) is minimum in the least-squares sense. A direct solution for u(v, τ′) requires computing the inverse of a large matrix which may have dimensions of 60 000 × 60 000 for a typical field data set (Section F.3).
To circumvent solving a problem that involves a very large matrix, Fourier transform the CMP gather d(h, t′) in the direction of the stretch variable t′. Correspondingly, apply Fourier transform to equation (13b) with respect to t′:
where ω′ is the Fourier dual of t′.
For each ω′, define d′ : d′(h, w′) and u : u(v, ω′) as complex vectors in h and v, respectively. Now, consider equation (14) in matrix notation:
The complex matrix L is given by equation (19) of Section F.3. The complex vectors d′ and u have lengths nh, the number of offsets, and nv, the number of constant velocities used in the transform defined by equation (14), respectively. The complex matrix L then has dimensions nh × nv. For a typical field data set, nh = 60 and nv = 60; hence, the complex matrix L may have dimensions of 60 × 60. As such, instead of solving one single problem using equation (13b) in the stretched time t′ domain that involves a very large matrix, we solve nω problems, where nω is the number of frequencies ω′, in the Fourier transform domain using equation (14) involving a small matrix L of equation (F-19).
We now restate our objective within the context of the matrix equation (15): For each ω′, estimate the complex vector u : u(v, ω′) such that the difference e : e(v, ω′) between the complex vector of the actual CMP gather d : d(h, ω′) and that of the modeled CMP gather d′ : d′(h, ω′) is minimum in the least-squares sense.
where T denotes transpose of the matrix L, the asterisk denotes complex conjugate and (LT*L)-1LT* is the least-squares (also called generalized) inverse of L.
Equation (16) gives the unconstrained least-squares solution for u, which may be unstable if the matrix LT*L does not have a stable inverse. Stability is attained by constraining the solution, which requires replacing LT*L in equation (16) with LT*L + βI, where β is called the damping factor (also called the Lagrange multiplier) and I is the identity matrix . Also, there can be singularities or near singularities in the matrix L. This is primarily because of the nonuniqueness of hyperbolic summation paths at and near zero offset, and discrete sampling over a finite range of offsets. The method of singular-value decomposition (SVD) is used to obtain the least-squares solution for u of equation (16) when the matrix L is singular or near singular (Section F.3).
We now outline the velocity-stack processing with reduced amplitude smearing based on the discrete hyperbolic Radon transform.
- Start with a CMP gather, d(h, t) and apply t2-stretching, d(h, t′ = t2).
- Fourier transform in the t′ direction, d(h, ω′).
- For each ω′, set up the L matrix (equation F-16) based on the geometry of the CMP gather and solve for u of equation (16) using the singular-value decomposition (Section F.3).
- Inverse Fourier transform to get u(v, τ′).
- Undo the t2-stretching to get u(v, τ), the velocity-stack gather with reduced amplitude smearing.
- Perform a desired operation, such as muting the zone of multiples, in the velocity-stack domain.
- Finally, perform inverse mapping back to the offset domain to get the modeled CMP gather d′(h, t) using equation (10b). During this inverse mapping, multiples, primaries, or all of the hyperbolic events can be modeled.
Recall that the minimum error associated with the least-squares solution is the difference between the actual CMP gather and the modeled CMP gather obtained from inverse mapping from the Radon transform domain back to the offset domain using equation (10b). Since the forward mapping defined by equation (10a) involves hyperbolic events, only, from the offset domain to the velocity domain, this difference gather should contain only the nonhyperbolic events, such as random or linear noise, that may be present in the original CMP gather.
- Beylkin, 1987, Beylkin, G., 1987, Discrete Radon transform: IEEE Trans. on Acoustics, Speech and Signal Processing, AASP-35, 162–172.
- Thorson and Claerbout (1985), Thorson, J. R. and Claerbout, J. F., 1985, Velocity-stack and slant-stack stochastic inversion: Geophysics, 50, 2727–2741.
- Yilmaz, 1989, Yilmaz, O., 1989, Velocity-stack processing: Geophys. Prosp., 37, 357–382.
- Lines and Treitel (1984), Lines, L. R. and Treitel, S., 1984, Tutorial: A review of least-squares inversion and its application to geophysical problems: Geophys. Prosp., 32, 159–186.
- Velocity-stack transformation
- The parabolic radon transform
- Practical considerations
- Impulse response of the velocity-stack operator
- Field data examples
- Radon-transform multiple attenuation