The discrete Radon transform

From SEG Wiki
Jump to navigation Jump to search
Seismic Data Analysis
Series Investigations in Geophysics
Author Öz Yilmaz
ISBN ISBN 978-1-56080-094-1
Store SEG Online Store

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 [1] 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.

Figure 6.4-2  (a) A synthetic CMP gather with three primary reflections; (b) a synthetic CMP gather with one primary reflection (arrival time at 0.2 s at zero-offset time) and its multiples; (c) composite CMP gather containing the primaries and multiples in (a) and (b); (d) the conventional velocity-stack gather derived from the composite CMP gather using equation (10a). Note the amplitude smearing along the velocity axis.

To reduce the amplitude smearing on conventional velocity-stack gathers Thorson and Claerbout [2] 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 [3]. Equation (9b) then takes the form


In the stretched coordinates, equations (10a, 10b) become:




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.

Following Lines and Treitel [4], the solution for equation (15) that minimizes the error vector e in the least-squares sense is derived in Section F.3 and is given by


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 [4]. 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.

  1. Start with a CMP gather, d(h, t) and apply t2-stretching, d(h, t′ = t2).
  2. Fourier transform in the t′ direction, d(h, ω′).
  3. 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).
  4. Inverse Fourier transform to get u(v, τ′).
  5. Undo the t2-stretching to get u(v, τ), the velocity-stack gather with reduced amplitude smearing.
  6. Perform a desired operation, such as muting the zone of multiples, in the velocity-stack domain.
  7. 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.

The velocity-stack gather estimated by equation (16) as described by the above sequence is one form of the generalized discrete forward Radon transform [1].

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.

Additional equations








  1. 1.0 1.1 Beylkin, 1987, Beylkin, G., 1987, Discrete Radon transform: IEEE Trans. on Acoustics, Speech and Signal Processing, AASP-35, 162–172.
  2. Thorson and Claerbout (1985), Thorson, J. R. and Claerbout, J. F., 1985, Velocity-stack and slant-stack stochastic inversion: Geophysics, 50, 2727–2741.
  3. Yilmaz, 1989, Yilmaz, O., 1989, Velocity-stack processing: Geophys. Prosp., 37, 357–382.
  4. 4.0 4.1 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.

See also

External links

find literature about
The discrete Radon transform
SEG button search.png Datapages button.png GeoScienceWorld button.png OnePetro button.png Schlumberger button.png Google button.png AGI button.png