# Multichannel filtering techniques for noise and multiple attenuation

Series Investigations in Geophysics Öz Yilmaz http://dx.doi.org/10.1190/1.9781560801580 ISBN 978-1-56080-094-1 SEG Online Store

### F.1 Analysis of guided waves

Marine data often are contaminated by guided waves that travel horizontally within the water layer or in the layers beneath the water layer. These waves exhibit characteristics that depend on water depth and on the geometry and material properties of the substrata. Modeling these pressure waves traveling within the water layer can lead to a better understanding of certain aspects of the field data and sometimes may even result in inferences about the strata below the water layer.

The well-known normal mode theory provides a way to laterally extrapolate acoustic and elastic waves [1] [2]. In this section, a normal mode procedure is applied to model shot profiles recorded over a water layer on top of a homogeneous elastic half space. Raypaths corresponding to multiple reflections, direct arrivals, refracted arrival, and its multiples are included in the normal mode theory.

The seismic waveguide effect of a surface layer is well known. Wave propagation in a surface layer, in particular guided waves, can be described by using the normal mode theory [1]. Pekeris’ model consists of a liquid layer over an acoustic (liquid) half space. More general models, which consist of a liquid layer on top of an elastic half space, were investigated by Press and Ewing [2]. The most complete summary of work in the field in the classic work by Ewing [3].

Guided waves are dispersive. This means each frequency component travels at a different speed; namely, the horizontal phase velocity. The dispersive character of guided waves is most pronounced in shallow water environments (less than 100 m). Depending on various water-bottom conditions, such as a mud layer with variable thickness or a hard bottom, the character of these waves may vary from shot to shot (Figure 6.0-3). They also can cause linear noise on stacked data (Figure 6.2-8a) and are easily confused with the linear noise that is associated with side scatterers (Figure 6.0-4).

McMechan and Yedlin [4] proposed a way to obtain phase velocity information from field data. Their approach is based on a wavefield transformation. The shot record first is transformed into the slant-stack domain. Fourier transforming (in time) each trace of the slant-stack gather then yields phase velocity as a function of frequency. This two-step process is demonstrated with the field data example in Figure F-1a. Line CC′ refers to the critical angle of propagation with a large reflection energy. To the right of this line is the supercritical propagation. The slant-stack gather is shown in Figure F-1b, and its 1-D amplitude spectrum is shown in Figure F-1c. Note that the horizontal axis in the slant-stack domain is the ray parameter, that is, the inverse of the horizontal phase velocity. Therefore, in Figure F-1c, we see the variation of the horizontal phase velocity as a function of frequency. Each curve corresponds to a particular normal mode propagating in the water layer. The phase velocities of the normal mode components asymptotically approach that of the water velocity vw at the high-frequency end of the spectrum.

Consider the recording geometry depicted in Figure F-2 with a source and a receiver cable, and subsurface model that consists of a water layer on top of an elastic half space. The source is at a certain depth below the water surface, so the two raypaths — primary and ghost, must be considered.

Figure F-1  (a) A shot gather containing the strong reflected and refracted multiples associated with hard water-bottom conditions. Here, CC′ = critical-angle energy. (b) The slant-stack gather derived from this shot gather. (c) The ω − p gather derived from the τ − p gather in panel (b). The inverse of p is the horizontal phase velocity. This figure demonstrates the dispersive nature of guided waves; that is, phase velocity is a function of frequency for all propagating normal mode components. These modes are represented by the curved trajectories on panel (c).

The normal mode theory of Ewing [3] provides an analytic expression for phase velocity as a function of frequency (the so-called characteristic equation or dispersion relation) for a given surface-layer model. The characteristic equation for the geometry in Figure F-2 is given by

 ${\displaystyle \tan(k_{x}Hr_{1})={\frac {\rho _{2}\beta _{2}^{4}r_{1}}{\rho _{1}c^{4}r_{2}}}\left[4r_{2}s_{2}-(1+s_{2}^{2})^{2}\right],}$ (F-1)

where kx is the horizontal wavenumber, which is equal to via equation (7), H is water depth, ρ1 and ρ2 are the water and substratum densities, respectively, β2 is the S-wave velocity of the substratum, and c is the phase velocity of the guided waves in the water layer. Finally, the normalized variables are

 ${\displaystyle k_{x}=p\omega .}$ (7)

 ${\displaystyle r_{1}={\sqrt {{\frac {c^{2}}{\alpha _{1}^{2}}}-1}},}$ (F-2a)

 ${\displaystyle r_{2}={\sqrt {{\frac {c^{2}}{\alpha _{2}^{2}}}-1}},}$ (F-2b)

 ${\displaystyle s_{2}={\sqrt {{\frac {c^{2}}{\beta _{2}^{2}}}-1}},}$ (F-2c)

where α1 is the P-wave velocity in the water layer and α2 is the P-wave velocity in the substratum. Because of the periodic nature of the tangent, equation (F-1) has a multiple-valued function on the left side and a single-valued function on the right side. To explicitly state the ambiguity, equation (F-1) can be rewritten as

 ${\displaystyle k_{x}Hr_{1}+n\pi =\tan ^{-1}B,}$ (F-3)

where B is the right-hand side of equation (F-1), and the integer n = 0, 1, 2, …, defines the mode number. By carefully examining equations (F-1) and (F-3), note that phase velocity c is a function of frequency ω; hence, guided waves are dispersive. Equation (F-3) yields real values of kx for α1cα2. (Here, we assume that α1 < β2 < α2.) Table F-1 gives a summary of the phase-velocity regions and types of rays associated with each region. Note that we get a wide range of wave types propagating each with a certain phase velocity or a range of phase velocities.

Figure F-2  The geometry for normal-mode modeling of the guided waves illustrated in Figure F-3. Here, S = the source, R = receivers, hs = source depth, and hr = receiver depth.

Only in the supercritical region, α1 < c < β2, are waves totally trapped within the water layer. These waves often form the major contribution to normal-mode propagation at longer offsets as in the field data example in Figure F-1a, where the supercritical region is to the right of line CC′. In the subcritical region, energy leaks into the substratum (thus, the name leaky modes). The contribution of this region to energy arriving at longer offsets is relatively weak.

The recorded pressure wavefield for the supercritical region at various receiver locations (Figure F-2) is given by

 ${\displaystyle P(x,z=h_{r},t)=4\sum _{n}\int d\omega \{\omega ^{2}A(\omega )\ \sin(k_{x}r_{1}h_{s})\ \sin(k_{x}r_{1}h_{r})\ \exp[i\omega (t-x\!/\!c)]\},}$ (F-4)

where A(ω) is the amplitude spectrum of the source. Yilmaz [5] modified this expression to account for the ghost effect. Note that reciprocity is satisfied here — the product of the sinusoidal factors that modulate the source spectrum A(ω) is unchanged if hr and hs are interchanged.

 Phase Velocity Ray Type α1 < c < α2 supercritical P-waves (wide-angle reflections) β1 < c < α2 supercritical, partially reflected P-waves c > α2 subcritical P-waves c = β2 critically refracted S-wave c = α2 critically refracted P-wave
Figure F-3  Superposition of all modes in water layers of different thicknesses. The depth model is shown in Figure F-2. Ghost effects are included in the modeling. Refer to the text for description of the labeled events.

We now consider modeling of normal modes using equation (F-4), including the ghost effect, for a range of water depths. The model parameters are: α1 = 1500 m/s, β2 = 2α1, α2 = 1.6β2, and ρ2/ρ1 = 2.2. All the experimental results represent impulse responses of guided waves — that is, A(ω) = 1 in equation (F-4). Guided waves in Figure F-3 manifest themselves with a complex interfering wave pattern at shallow water, then gradually separate into simple water-bottom multiples at increasing water depths. The dispersive character of the guided waves is prominent, especially for shallow water depths. In Figure F-3, the guided waves are simulated in the supercritical region. The elastic substratum, which is equivalent to the hard water-bottom case, supports the early refraction energy RP, its multiple RM, and the reflected water-bottom multiples M1, M2, and M3. The acoustic substratum (β2 = 0, equivalent to the soft water-bottom case), yields only the reflected water-bottom multiples. Acoustic behavior of the substratum implies that no P-to-S conversion occurs. The phase velocity curves in Figure F-4 verify the existence of a number of propagating modes for each case. Note the cable truncation effect (CT).

### F.2 Wavefield extrapolation in the τ − p domain

In the slant-stack transform, we discussed slant-stack transformation of a wavefield from midpoint-offset to midpoint-ray-parameter coordinates. This transformation is done by applying linear moveout and summing over the offset range for each value of the ray-paramater. From the results in Section D.1, the double square-root operator can be specialized for migration before stack in midpoint-ray-parameter coordinates [6]. To derive the extrapolation equation in the slant-stack domain, we start with the relationships between the transform-domain variables:

 ${\displaystyle p={\frac {\sin \theta }{v}}}$ (F-5a)

and

 ${\displaystyle {\frac {vk_{h}}{2\omega }}=\sin \theta ,}$ (F-5b)

where kh is the offset wavenumber, ω is the temporal frequency, p is the ray parameter, v is the propagation velocity, and θ is the angle of propagation measured from the vertical. The normalized offset wavenumber H is defined as (Section D.1)

 ${\displaystyle H={\frac {vk_{h}}{2\omega }}.}$ (F-5c)

Combining the relationships given by equations (F-5a,F-5b) with the definition given by equation (F-5c), we get

 ${\displaystyle H=pv.}$ (F-6)
Figure F-4  Phase velocity as a function of frequency for the cases shown in Figure F-3. Here, CT denotes the cable truncation artifact.

We rewrite the double square-root equation from Section D.1

 ${\displaystyle DSR(Y,H)={\sqrt {1-(Y+H)^{2}}}+{\sqrt {1-(Y-H)^{2}}},}$ (F-7)

and substitute the definition of H in terms of p from equation (F-6) to get

 ${\displaystyle DSR(Y,pv)={\sqrt {1-(Y+pv)^{2}}}+{\sqrt {1-(Y-pv)^{2}}}.}$ (F-8)

Ottolini [6] used this operator for migration before stack in midpoint-ray-parameter coordinates. The procedure is described in Figure F-5.

Clayton and McMechan [7] adapted equation (F-8) to zero-dip case which is equivalent to setting Y = 0:

 ${\displaystyle DSR(Y=0,H=pv)=2{\sqrt {1-p^{2}v^{2}}}.}$ (F-9)

They then used this operator to downward continue refracted waves on CMP or shot gathers. The objective of inversion of a refraction profile is to estimate a velocity profile in depth, v(z). The procedure is outlined in Figure F-6. The final step yields a profile, P(p, z), of horizontal phase velocity 1/p as a function of depth z. Two issues must be kept in mind. First, the procedure is based on a layered earth assumption. Second, the procedure requires knowledge of the medium velocity to extrapolate the wavefield in depth. To get around this second problem, the process must be iterated until the phase velocity profile converges to the velocity function used in the extrapolation. It usually takes up to three iterations to achieve convergence.

Figure F-5  Flowchart for migration before stack in midpoint-ray-parameter coordinates.

### F.3 Mathematical foundation of the discrete radon transform

The forward Radon transform u(r, τ) of a 2-D continuous function d(h, t) is given by the integral expression [8]

 ${\displaystyle u(r,\tau )=\int _{-\infty }^{\infty }d\,[h,t=\tau +\phi (r,h)]\ dh,}$ (F-10a)

where h and t are the input variables, and r and τ are the transform variables. The integration is along trajectories expressed as linear functions of traveltimes t and τ.

The inverse Radon transform given by the integral expression

 ${\displaystyle d(h,t)=\int _{-\infty }^{\infty }\rho (\tau )\ast u[r,\tau =t-\phi (r,h)]\ dr}$ (F-10b)

incorporates convolution of u(r, τ) with the rho filter ρ(τ) prior to integration. For the 2-D data type d(h, t), the Fourier transform of the rho filter is of the form ${\displaystyle {\sqrt {\omega }}\ \exp(i\pi \!/\!4).}$

In practice, of course, we do not deal with continuous functions; instead, we have discretely sampled data in time and space. So we need to replace the integrals with discrete summations in equations (F-10a,F-10a). We can use the least-squares technique to compute the discrete Radon transform [8] and account for the effects of discrete sampling and finite cable length. We shall consider three different forms of the function ϕ(r, h) in equation (F-10a) and, accordingly, define discrete hyperbolic [9] [10]. [11], parabolic [12], and linear [9] Radon transforms.

Figure F-6  Processing flow for inversion of refraction profiles.

First, we shall discuss the discrete hyperbolic Radon transform. The transform variable r of the function ϕ(r, h) in equation (F-10a) represents stacking velocity v. The discrete transformation of the CMP data from the offset domain to the velocity domain is achieved by hyperbolic moveout correction and summing over offset

 ${\displaystyle u(v,\tau )=\sum _{h}d(h,t={\sqrt {\tau ^{2}+4h^{2}\!/\!v^{2}}}),}$ (F-11a)

where t is the two-way traveltime, τ is the two-way zero-offset time, and h is the half-offset. The inverse transform from the velocity domain back to the offset domain is achieved by inverse hyperbolic moveout correction and summing over velocity

 ${\displaystyle d'(h,t)=\sum _{v}u(v,\tau ={\sqrt {t^{2}-4h^{2}\!/\!v^{2}}}).}$ (F-11b)

The relationship between (h, t) and (v, τ) coordinates is given by the hyperbolic moveout equation:

 ${\displaystyle t^{2}=\tau ^{2}+{\frac {4h^{2}}{v^{2}}}.}$ (F-12a)

To comply with the linear form of the Radon transform defined by equation (F-11a), apply stretching in the time direction by setting t′ = t2 and τ′ = τ2. Equation (F-12a) then takes the form

 ${\displaystyle t'=\tau '+{\frac {4h^{2}}{v^{2}}}.}$ (F-12b)

Also in the stretched coordinates, equations (F-11a, F-11b) take the forms

 ${\displaystyle u(v,\tau ')=\sum _{h}d(h,t'=\tau '+4h^{2}\!/\!v^{2})}$ (F-13a)

and

 ${\displaystyle d'(h,t')=\sum _{v}u(v,\tau '=t'-4h^{2}\!/\!v^{2}).}$ (F-13b)

Now, consider equation (F-13b) in matrix notation:

 ${\displaystyle \mathbf {d} '=\mathbf {Lu} .}$ (F-14)

The matrix operator L maps each point in u : u(v, τ′) onto a parabola in d′ : d′(h, t′), the modeled CMP gather in the stretched coordinates.

The objective is to estimate a u : u(v, τ′) such that the difference e : (h, t′) between the actual CMP gather d : d(h, t′) and the modeled CMP gather d′ : d′(h, t′) is minimum in the least-squares sense. Using the matrix notation and equation (F-14), e is defined as

 ${\displaystyle \mathbf {e} =\mathbf {d} -\mathbf {Lu} .}$ (F-15)

To distinguish from the conventional velocity-stack gather defined by equation (F-11a) with amplitude smearing, we shall refer to u in equation (F-14) as the Radon transform of d in the stretched coordinates. By the Radon transform, parabolas in the offset domain (h, t′) with the stretched coordinate t′ (or, equivalently, hyperbolas in the offset domain h, t) with the unstretched coordinate t are represented by points in the velocity domain (v, τ).

Following Lines and Treitel [13], the least-squares solution for equation (F-14) is determined, first, by expressing the cumulative squared error S as

 ${\displaystyle S=\mathbf {e^{T}e} .}$ (F-16a)

where T is for transpose. By substituting for e from equation (F-15), we get

 ${\displaystyle S=(\mathbf {d} -\mathbf {Lu} )^{\mathbf {T} }(\mathbf {d} -\mathbf {Lu} ).}$ (F-16b)

Minimization of S with respect to u yields the desired least-squares solution:

 ${\displaystyle \mathbf {u} =(\mathbf {L^{T}L} )^{\mathbf {-1} }\mathbf {L^{T}d} ,}$ (F-17)

where T denotes transpose of the matrix L, and (LTL)-1LT is the least-squares (also called generalized) inverse of L.

In equation (F-17) d is a column vector containing all the data points from the actual CMP gather in the stretched coordinates. Its length is nhnt, where nh is the number of offsets and nt is the number of time samples in t′. Also, u is a column vector containing all the points from the velocity stack-gather; its length is nvnτ, where nv is the number of velocities and nτ is the number of time samples in τ′. The elements of the matrix operator L are delta functions in the velocity domain; its dimensions are nhnt × nvnτ. For a typical field data set, nh = 60, nt = 1000, nv = 60 and nτ = 1000; this implies an L matrix of dimension 60 000 × 60 000.

The minimum error e of equation F-15 associated with the least-squares solution u of equation (F-17) should be interpreted as being the CMP gather that contains only the nonhyperbolic events, such as random or linear noise, that may be present in the original CMP gather d. The velocity-stack gather u obtained from equation (F-17) is one form of the discrete Radon transform of d [8].

Direct solution for u of equation (F-17) requires computing the inverse of the matrix LTL which may have dimensions of 60 000 × 60 000 for a typical field data set. Inverting such a large matrix is quite impractical. To circumvent solving a problem that involves a very large matrix L, Fourier transform the CMP gather d(h, t′) in the direction of the stretch variable t′. Correspondingly, apply Fourier transform to equation (F-13b) with respect to t′ given by

 ${\displaystyle d'(h,\omega ')=\sum _{v}u(v,\omega ')\ \exp(-i\omega '4h^{2}\!/\!v^{2}),}$ (F-18)

where ω′ is the Fourier dual of t′. For each ω′ component of d′(h, ω′) and u(v, ω′), equation (F-18) can be written in the matrix form of equation (F-15), where L now is a complex matrix of the form

 ${\displaystyle \mathbf {L} ={\begin{pmatrix}e^{-i\omega '4h_{1}^{2}/v_{1}^{2}}&e^{-i\omega '4h_{1}^{2}/v_{2}^{2}}&\ldots &e^{-i\omega '4h_{1}^{2}/v_{n}^{2}}\\e^{-i\omega '4h_{2}^{2}/v_{1}^{2}}&e^{-i\omega '4h_{2}^{2}/v_{2}^{2}}&\ldots &e^{-i\omega '4h_{2}^{2}/v_{n}^{2}}\\{\boldsymbol {\vdots }}&{\boldsymbol {\vdots }}&{\boldsymbol {\ddots }}&{\boldsymbol {\vdots }}\\e^{-i\omega '4h_{m}^{2}/v_{1}^{2}}&e^{-i\omega '4h_{m}^{2}/v_{2}^{2}}&\ldots &e^{-i\omega '4h_{m}^{2}/v_{n}^{2}}\end{pmatrix}},}$ (F-19)

with dimensions m × n = nh × nv, where nh is the number of offsets, nv is the number of velocities, and d′ and u are complex vectors of lengths nh and nv, respectively. Note that the elements of the L matrix depend on the geometry of the input CMP gather and the range of velocities used in constructing the velocity-stack gather.

We now restate our objective in solving for u of the matrix equation (F-14) within the context of the Fourier transform domain: 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.

As for equation (F-17), follow the steps involving equations (F-14) through (F-16) to derive the least-squares solution in the Fourier transform domain as

 ${\displaystyle \mathbf {u} =(\mathbf {L^{T\ast }L} )^{\mathbf {-1} }\mathbf {L^{T\ast }d} ,}$ (F-20)

where the asterisk denotes complex conjugate.

Equation (F-20) gives the unconstrained least-squares solution for u. To avoid singularities or near singularities in the matrix LT*L, the solution is constrained by incorporating a damping factor β (also called the Lagrange multiplier) into equation (F-20) [13]:

 ${\displaystyle \mathbf {u} =(\mathbf {L^{T\ast }L} +\beta \mathbf {I} )^{\mathbf {-1} }\mathbf {L^{T\ast }d} .}$ (F-21)

Because of the near-singular character of the complex matrix L, especially for small values of ω′, the solution given by equation (F-21) is best reformulated in terms of the singular-value decomposition (SVD) of the matrix L [14]. This procedure factors L into a product of three matrices:

 ${\displaystyle \mathbf {L} =\mathbf {U\Lambda V^{T\ast }} .}$ (F-22)

By using this factorized form of the matrix L, the constrained solution given by equation (F-21) takes the form

 ${\displaystyle \mathbf {u} =\mathbf {V} [(\mathbf {\Lambda ^{2}} +\beta \mathbf {I} )^{\mathbf {-1} }\mathbf {\Lambda } ]\mathbf {U^{T\ast }d} ,}$ (F-23)

where

 ${\displaystyle (\mathbf {\Lambda ^{2}} +\beta \mathbf {I} )^{\mathbf {-1} }\mathbf {\Lambda } ={\begin{pmatrix}\Gamma _{1}&0&\ldots &0\\0&\Gamma _{2}&\ldots &0\\\vdots &\vdots &\ddots &\vdots \\0&0&\ldots &\Gamma _{n_{v}}\end{pmatrix}},}$ (F-24)

and ${\displaystyle \Gamma _{i}=\lambda _{i}\!/\!(\lambda _{i}^{2}+\beta ),}$ with λi as the positive square roots of the eigenvalues ${\displaystyle \lambda _{i}^{2}}$ of LTL. Recall that the damping factor β is a scalar that prevents the solution (equation F-23) from becoming unstable.

We now summarize the method for computing the Radon transform represented by a velocity-stack gather:

1. Start with a CMP gather d(h, t) and apply the t2-stretching, d(h, t′ = t2).
2. Fourier transform in the t′ direction, d(h, ω′).
3. For a specific value of ω′:
1. Set up the L matrix (equation F-19) based on the geometry of the CMP gather.
2. Set up the d vector by transposing the data set d(h, ω′).
3. Apply singular-value decomposition to L equation (F-22), and compute U, Λ and VT*, hence UT* and V.
4. Specify a value for the damping factor β and set up the diagonal matrix of equation (F-24).
5. Finally, solve for u (equation F-23).
6. Repeat (3) for all ω′ values and accumulate the results in u(v, ω′).
1. Inverse Fourier transform to get u(v, τ′).
2. Undo the t2-stretching to get u(v, τ); this is the desired result, namely the velocity-stack gather with reduced amplitude smearing, which is a special form of the discrete Radon transform of d(h, t)
3. The inverse mapping back to the offset domain to get the modeled CMP gather is done using equation (F-11b).

An important practical question is the sampling along the velocity axis in the transform domain. Specifically, in equation (F-19), one may consider sampling in terms of Δv or 1/Δv2. This will affect the quality of the velocity stack gather and the reconstructed gather. When 1/Δv2 sampling is used, the structure of the L matrix in equation (F-19) becomes Toeplitz [15]. Hence, one can use Levinson recursion to solve equation (F-20) and thus obtain the Radon transform u.

Figure F-7 shows a comparison of the performance of the various solutions to equation (F-20). The input is a synthetic CMP gather with band-limited noise. The velocity-stack gather was computed using the SVD approach described here with Δv and 1/Δv2 sampling criteria and the Toeplitz solution with the 1/Δv2 sampling criterion. Figure F-7 also shows the reconstructed gathers using equation (F-11b). Results indicate that the SVD approach with the Δv sampling criterion best preserves the amplitudes at all offsets. The 1/Δv2 sampling causes loss of amplitudes at near offsets.

Figure F-7  Solutions to equation (F-20) using SVD and Toeplitz procedures with two different sampling criteria — Δv and 1/Δv2. (a) A synthetic CMP gather as in Figure 6.4-2c with band-limited noise; (b), (c) and (d) the velocity-stack gathers using three combinations of sampling and procedure as indicated on top of each panel; and (e), (f) and (g) reconstructed gathers from (b), (c) and (d), respectively.

We now set up the least-squares problem to perform the discrete parabolic Radon transform [12], [16]. As an alternative to stretching in the time direction as defined by equation (F-12b) to comply with the linear form of the Radon transform (equation F-11a), consider the application of normal moveout correction to input CMP data (equation F-12a) using a velocity function vn:

 ${\displaystyle t_{n}={\sqrt {t^{2}-{\frac {4h^{2}}{v_{n}^{2}}}}},}$ (F-25)

such that hyperbolic events in the original gather d(h, t) can be approximately represented by parabolic events with the traveltime equation

 ${\displaystyle t_{n}=\tau +qh^{2},}$ (F-26)

where tn is the time after NMO correction, vn is the hyperbolic moveout correction velocity function, τ is the two-way zero-offset time, and q is the parameter that defines the curvature of the parabolic trajectory described by equation (F-26).

In the coordinates of the NMO-corrected gather d(h, tn), equations (F-11a, F-11b) take the forms

 ${\displaystyle u(q,\tau )=\sum _{h}d(h,t_{n}=\tau +qh^{2})}$ (F-27a)

and

 ${\displaystyle d'(h,t_{n})=\sum _{q}u(q,\tau =t_{n}-qh^{2}).}$ (F-27b)

We want to find an estimate of u such that when inverse transformed, the difference between the modeled moveout-corrected CMP gather d′(h, tn) and the original moveout-corrected CMP gather d(h, tn), which is the desired output, is minimum in the least-squares sense.

A practical estimation of u based on the least-squares scheme can be achieved, again, by first Fourier transforming equation (F-26b) as

 ${\displaystyle d'(h,\omega _{n})=\sum _{q}u(q,\omega _{n})\ \ \exp(-i\omega _{n}qh^{2}),}$ (F-28)

where ωn is the Fourier dual of tn. For each ωn component of d′(h, ωn) and u(q, ωn), equation (F-28) can be written in the matrix form of equation (F-14), where L is a complex matrix of the form:

 ${\displaystyle \mathbf {L} ={\begin{pmatrix}e^{-i\omega _{n}q_{1}h_{1}^{2}}&e^{-i\omega _{n}q_{2}h_{1}^{2}}&\ldots &e^{-i\omega _{n}q_{n}h_{1}^{2}}\\e^{-i\omega _{n}q_{1}h_{2}^{2}}&e^{-i\omega _{n}q_{2}h_{2}^{2}}&\ldots &e^{-i\omega _{n}q_{n}h_{2}^{2}}\\{\boldsymbol {\vdots }}&{\boldsymbol {\vdots }}&{\boldsymbol {\ddots }}&{\boldsymbol {\vdots }}\\e^{-i\omega _{n}q_{1}h_{m}^{2}}&e^{-i\omega _{n}q_{2}h_{m}^{2}}&\ldots &e^{-i\omega _{n}q_{n}h_{m}^{2}}\\\end{pmatrix}},}$ (F-29)

with dimensions m × n = nh × nq, where nh is the number of offsets, and nq is the number of q parameters for which the Radon transform u(q, τ) is to be estimated. As for the hyperbolic Radon transform (equation F-19), note that the elements of the L matrix depend only on the geometry of the input CMP gather.

Aside from the L matrix, for each frequency component ωn, define the vectors for the input gather d : d(h, t), the unknown transform u : u(q, τ), and the modeled CMP gather d′ : d′(h, t). As for equation (F-18), for each frequency component ωn, equation (F-28) is written in the matrix form of equation (F-14), where L is the complex matrix given by equation (F-29). Using this new form of the L matrix, follow the same procedure starting with equation (F-20) to obtain the solution expressed by equation (F-23).

Finally, we set up the least-squares problem to perform the discrete linear Radon transform. The transform variable r of the function ϕ(r, h) in equation (F-11a) represents the ray parameter p. The relationship between (h, t) and (p, τ) coordinates is given by the parabolic moveout equation

 ${\displaystyle t=\tau +2ph.}$ (F-30)

The transformation of the CMP data from the offset domain to the ray-parameter domain is achieved by applying linear moveout correction and summing over offset

 ${\displaystyle u(p,\tau )=\sum _{h}d(h,t=\tau +2ph),}$ (F-31a)

where t is the two-way traveltime, τ is the intercept time at p = 0, and h is the half-offset. The inverse transform from the ray-parameter domain back to the offset domain is achieved by inverse linear moveout correction and summing over velocity

 ${\displaystyle d'(h,t)=\sum _{p}u(v,\tau =t-2ph).}$ (F-31b)

We want to estimate u : u(p, τ) such that when inverse transformed back to the offset domain, the difference between the modeled CMP gather d′ : d′(h, t) and the original CMP gather d : d(h, t) is minimum in the least-squares sense.

A practical estimation of u based on the least-squares scheme can be achieved, once again, by first Fourier transforming equation (F-31b) as

 ${\displaystyle d'(h,\omega )=\sum _{p}u(p,\omega )\ \ \exp(-i\omega 2ph).}$ (F-32)

For each ω component of d′(h, ω) and u(p, ω), equation (F-32) can be written in the matrix form of equation (F-14), where L is a complex matrix of the form

 ${\displaystyle \mathbf {L} ={\begin{pmatrix}e^{-i\omega 2p_{1}h_{1}}&e^{-i\omega 2p_{2}h_{1}}&\ldots &e^{-i\omega 2p_{n}h_{1}}\\e^{-i\omega 2p_{1}h_{2}}&e^{-i\omega 2p_{2}h_{2}}&\ldots &e^{-i\omega 2p_{n}h_{2}}\\{\boldsymbol {\vdots }}&{\boldsymbol {\vdots }}&{\boldsymbol {\ddots }}&{\boldsymbol {\vdots }}\\e^{-i\omega 2p_{1}h_{m}}&e^{-i\omega 2p_{2}h_{m}}&\ldots &e^{-i\omega 2p_{n}h_{m}}\\\end{pmatrix}},}$ (F-33)

with dimensions m × n = nh × np, where nh is the number of offsets, and np is the number of ray parameters p for which the Radon transform u(p, τ) is to be estimated. As for the hyperbolic and parabolic cases, the elements of the L matrix depend only on the geometry of the input CMP gather.

Aside from the L matrix, for each frequency component ω, define the vectors for the input gather d : d(h, t), the unknown transform u : u(p, τ), and the modeled CMP gather d′ : d′(h, t). As for the hyperbolic and parabolic cases, for each frequency component ω, equation (F-32) is written in the matrix form of equation (F-14), where L is the complex matrix given by equation (F-33). Using this new form of the L matrix, follow the same procedure starting with equation (F-20) to obtain the solution expressed by equation (F-23).

Figure F-8 shows a CMP gather and the modeled CMP gathers computed by inverse transforming the velocity-stack and slant-stack gathers that were themselves estimated using the respective discrete Radon transforms. The reconstruction based on the velocity-stack gather restores that component of the original gather associated with the hyperbolic events. The reconstruction based on the slant-stack transformation, when implemented as a special form of the discrete Radon transform, restores virtually all of the data characteristics of the original gather.

There exists an affinity between slant-stack transformation (the slant-stack transform) and the discrete linear Radon transform discussed here. To compute the slant-stack gather, recall that we applied linear moveout to input data and summed over the offset axis (equation 5). That operation is equivalent to applying the LT* matrix on the input vector d in equation (F-23). We then applied the rho filter to compensate for the attenuation of high frequencies during the summation over the offset axis. Applying the rho filter is to some extent equivalent to applying the additional matrix operator (LT*L)-1 of equation (F-20). Actually this 2-D operator does more than what the 1-D rho filter does — it accounts for the discrete sampling along the spatial axis and the finite cable length. Nevertheless, in practice, it turns out that just applying the rho filter in lieu of applying this computationally more involved operator yields a reasonably accurate slant-stack gather.

Figure F-8  (a) A CMP gather; (b) reconstructed CMP gather by way of velocity-stack transformation; (c) reconstructed CMP gather by way of slant-stack transformation. Both velocity-stack and slant-stack transforms are special forms of the discrete Radon transform.

### F.4 Free-surface multiple attenuation

Refer to the raypath configurations for the various types of multiples depicted in Figure 6.0-18. Most significant multiples involve one or more bounces from the free surface. Figure 6.0-18a shows the first- and second-order water-bottom multiples, and Figure 6.0-18b shows the first- and second-order free-surface multiples associated with a deeper reflector. Figure 6.0-18c shows the peg-leg multiples associated with a primary reflection and intrabed multiple reflection, and Figure 6.0-18d shows the first- and second-order intrabed multiple reflections. Finally, Figure 6.0-18e shows the first- and second-order interbed multiple reflections.

Figure F-9  (a) Marine recording geometry, (b) Noah’s geometry. S is the vertically downgoing source wavefield, X is the recorded upcoming wavefield with and Y is the wavefield without free-surface multiples. All wavefields are assumed to be plane waves. Adapted from Riley and Claerbout [17].

Shown in Figure 6.0-18 are but a limited subset of a countless group of multiples that are normally present in recorded marine data. Rather than attempting to attenuate all multiples, it makes sense to develop a technique aimed at attenuating the most significant set of multiples — those associated with the free surface. Riley and Claerbout [17] developed an elegant theory for free-surface multiple attenuation applicable to 1-D seismograms. They also extended their ideas to 2-D recorded data to attenuate diffracted multiples. The free-surface multiple attenuation theory by Riley and Claerbout [17] was reworked and widely publicized later by Kennett [18], Verschuur [19], Verschuur [20], Carvalho [21], Dragoset and McKay [22], Verschuur and Berkhout [23], and Berkhout and Verschuur [24]. Applications to field data were reported by Dragoset and McKay [22], Verschuur [25], Kelamis and Verschuur [26], and Dragoset and Jericevic [27]. Here, we shall develop the theory of free-surface multiple attenuation based on the work by Riley and Claerbout [17], then review the work by Verschuur [20].

Consider the 1-D earth model in Figure F-9a and recording using a source represented by a vertically downward traveling plane wave W. If the free surface were removed from the earth model (Figure F-9b), then the associated multiples would also be removed from the recorded 1-D seismogram. The hypothetical recording geometry that facilitates removal of the free surface is called Noah’s geometry by Riley and Claerbout [17], and the technique to remove the free-surface multiples based on Noah’s geometry is called Noah’s deconvolution. Note that in Noah’s geometry the recording cable now is being towed by a submarine rather than a surface vessel since the water level has risen as a result of the flood.

Whether you record with Noah’s geometry (Figure F-9b) or the actual geometry at the same datum level (Figure F-9a), the earth’s response E(z) to the same source S(z) defined in the z-transform domain should be the same. The earth’s response E(z) can be described by the ratio of the upcoming waves U(z) to the downgoing waves D(z):

 ${\displaystyle E(z)={\frac {U(z)}{D(z)}}.}$ (F-34a)

From Figure (F-9a), the upcoming wave is U(z) = −X(z), the recorded seismogram, and the downgoing wave is D(z) = S(z) + X(z). Here, we assume that the reflection coefficient of the free surface is −1. From Figure F-9b, the upcoming wave is U(z) = −Y(z), the seismogram that is free of surface multiples, and the downgoing wave is D(z) = S(z). Hence, the earth’s response E(z) defined by equation (F-34a) for the geometry of Figure F-9a is

 ${\displaystyle E(z)=-{\frac {X(z)}{S(z)+X(z)}},}$ (F-34b)

and for the geometry of Figure (F-9b) is

 ${\displaystyle E(z)=-{\frac {Y(z)}{S(z)}}.}$ (F-34c)

Combine equations (F-34b, F-34c) to obtain

 ${\displaystyle {\frac {Y(z)}{S(z)}}={\frac {X(z)}{S(z)+X(z)}}.}$ (F-35)

Assume that the source waveform is minimum phase (inverse filtering) and define its inverse F(z) = S−1(z). Then, solve equation (F-35) for the seismogram Y(z) free of surface multiples as

 ${\displaystyle Y(z)={\frac {X(z)}{1+F(z)X(z)}}.}$ (F-36)

Equation (F-38) is the basis for Noah’s deconvolution described by Riley and Claerbout [17]. A practical perspective on the use of equation (F-36) is gained by the series expansion of the right-hand side:

 ${\displaystyle Y(z)=X(z)[1-F(z)X(z)+F^{2}(z)X^{2}(z)-\cdots ],}$ (F-37)

where it is assumed that F(z)X(z) ≪ 1.

Equation (F-37) can be put into the convolutional form

 ${\displaystyle y(t)=x(t)-x(t)\ast [f(t)\ast x(t)]+x(t)\ast [f(t)\ast x(t)]\ast [f(t)\ast x(t)]-\cdots .}$ (F-38)

An iterative scheme to estimate the seismogram y(t) free of surface multiples can be devised based on equation (F-38). First, write the following recursive relation:

 ${\displaystyle y_{i+1}(t)=x(t)-x(t)\ast [f(t)\ast y_{i}(t)],\qquad i=0,1,2,\ldots .}$ (F-39)

where y0(t) = x(t).

The recursive relation described by equation (F-39) yields the truncated form of the series on the right-hand side of equation (F-38). To see this, set i = 0 in equation (F-39):

 ${\displaystyle y_{1}(t)=x(t)-x(t)\ast [f(t)\ast x(t)],}$ (F-40a)

which corresponds to the first two terms of the series in equation (F-39). Next, set i = 1 in equation (F-39) and substitute equation (F-40a) for y1 to get

 ${\displaystyle y_{2}(t)=x(t)-x(t)\ast [f(t)\ast x(t)]+x(t)\ast [f(t)\ast x(t)]\ast [f(t)\ast x(t)],}$ (F-40b)

which corresponds to the first three terms in equation (F-38).

Thus, by recursive application of equation (F-39) for a specified number of iterations, an estimate of the seismogram y(t) ≈ yi+1(t) that is free of most of the surface-related multiples can be obtained. Specifically, equation (F-39) states that the output yi+1(t) of the next iteration is the difference between the input trace x(t) and the input trace convolved with the inverse of the source wavelet f(t) convolved with the output of the previous iteration yi(t). The number of iterations in practice often is no more than three.

While the recursive scheme described here based on the work by Riley and Claerbout [17] is designed for a zero-offset seismogram x(t) associated with a 1-D earth model, the theory has been extended to handle nonzero-offset data associated with a 2-D earth model by the same authors, and subsequently by Verschuur [19], Verschuur and Berkhout [23], and Dragoset and McKay [22]. In fact, the recursive relation given by equation (F-39), with a group of shot records as input, is the same as the update formula given by Verschuur and Berkhout [23].

To attenuate free-surface multiples caused by reflectors with complex geometry, in theory, you need all of the shot records along the seismic line as one single input to equation (F-39) [22]. In practice, however, assumptions are made about the subsurface complexity to be able to use individual CMP gathers as input [26]. Also, assuming that spiking deconvolution has been applied to the input CMP gather prior to free-surface multiple attenuation, the inverse of the source wavelet f(t) in equation (F-39) may be set to a unit spike at t = 0. Then, it follows that

 ${\displaystyle P_{i+1}(h,t)=P_{0}(h,t)-P_{0}(h,t)\ast P_{i}(h,t),}$ (F-41)

where P0(h, t) is the input CMP gather and Pi(h, t) is the gather after free-surface multiple attenuation, and h and t are the offset and time variables, respectively.

The recursion described by equation (F-41) requires 2-D convolution of the input gather P0(h, t) with the result of the previous iteration Pi(h, t). For computational efficiency, the recursion can be performed in the space-frequency (h − ω) domain

 ${\displaystyle P_{i+1}(h,\omega )=P_{0}(h,\omega )-P_{0}(h,\omega )P_{i}(h,\omega ),}$ (F-42)

where P0(h, w) is the input CMP gather and Pi(h, ω) is the gather after free-surface multiple attenuation, both in the h − ω domain.

In practice, there are two issues that hinder the effectiveness of the free-surface multiple attenuation technique discussed here. First is the problem of an accurate estimate of the source wavelet s(t). Actually, it is the inverse of the source wavelet f(t) that is needed in the recursive relation described by equation (F-39). Deterministically, the source wavelet may be assumed equivalent to the recorded far-field signature (field data examples). Statistically, it may be assumed to be the minimum-phase inverse of the spiking deconvolution operator (inverse filtering). It is this latter viewpoint that is often taken in practice. As such, the input to the recursive estimate of the surface-related multiples in equation (F-41) is the CMP gather with spiking deconvolution.

The second issue is the missing near offsets. Attenuation of free-surface multiples using the recursive relation of equation (F-41) involves predicting multiples from the primaries. The free-surface multiples generated by a primary reflection that occurs within the near-offset range that is not recorded then are not attenuated. In fact, the problem of missing near-offset data has the most significant detrimental impact on the ability of the technique to attenuate all possible surface-related multiples [27]. The Radon transform provides a theoretically appealing way to estimate the unrecorded near offsets (the radon transform).

Figure F-10 shows a field data example of free-surface multiple attenuation. The free-surface multiples in the shot record (Figure F-10a) include the water-bottom multiples and the peg-leg multiples associated with the top-salt and base-salt reflections. The same record after free-surface multiple attenuation (Figure F-10b) exhibits some residual of these multiples, primarily because of the practical issues discussed above. The free-surface multiples heavily contaminate the subsalt reflections (Figure F-10c). After multiple attenuation, the subsalt region is largely free of surface-related multiples.

### F.5 Water-bottom multiple attenuation

A subset of surface-related multiples is water-bottom multiples. Based on the work by Morley [28], Wiggins [29] formulated a technique to attenuate water-bottom multiples associated with a complex water-bottom geometry. As for the free-surface multiple attenuation technique, this method also predicts multiples from the primaries contained in the recorded data within the cable length. To develop the conceptual basis, here, we shall review the technique for the special case of flat water bottom.

Refer to the raypath configuration in Figure F-11. Assume a horizontally layered earth model and treat a CMP gather a single wavefield equivalent to a shot record. The downgoing wavefield D(h, zW, ω) at the water bottom is given by

 ${\displaystyle D(h,z_{W},\omega )=P_{0}(h,z_{R},\omega )\ \exp[ik_{z}(z_{W}+z_{R})],}$ (F-43)

where P0(h, zR, ω) is the CMP gather that represents the recorded wavefield at the receiver depth zR, h is the offset variable, zW is the water depth, ω is the temporal frequency and kz is the vertical wavenumber. The upcoming wavefield U(h, zW, ω) at the water bottom is given by

 ${\displaystyle U(h,z_{W},\omega )=P_{0}(h,z_{R},\omega )\ \exp[-ik_{z}(z_{W}-z_{R})].}$ (F-44)

In equations (F-43) and (F-44), the sign convention defined in Section D.1 for wavefield extrapolation is used.

The CMP gather Pp(h, zW, ω) free of water-bottom multiples defined at the water bottom is given by [29]

 ${\displaystyle P_{p}(h,z_{W},\omega )=U(h,z_{W},\omega )-r_{W}(\omega )D(h,z_{W},\omega ),}$ (F-45)

where rW(ω) is the water-bottom reflectivity.

Substitute equations (F-43) and (F-44) into equation (F-45) and simplify to get

 ${\displaystyle P_{p}(h,z_{W},\omega )=P_{0}(h,z_{R},\omega )\{\ \exp[-ik_{z}(z_{W}-z_{R})]-r_{W}(\omega )\ \exp[ik_{z}(z_{W}+z_{R})]\}.}$ (F-46)

The CMP gather Pp(h, zR, ω) free of water-bottom multiples defined at the receiver depth is computed by extrapolating the CMP gather Pp(h, zW, ω) defined at the water bottom:

 ${\displaystyle P_{p}(h,z_{R},\omega )=P_{p}(h,z_{W},\omega )\ \exp[ik_{z}(z_{W}-z_{R})].}$ (F-47)

Now substitute equation (F-46) into (F-47) and simplify to get

 ${\displaystyle P_{p}(h,z_{R},\omega )=P_{0}(h,z_{R},\omega )[1-r_{W}(\omega )\ \exp(ik_{z}2z_{W})].}$ (F-48)

Finally, rewrite this equation as

 ${\displaystyle P_{p}(h,\omega )=P_{0}(h,\omega )-r_{W}(\omega )P_{0}(h,\omega )\ \exp(ik_{z}2z_{W}),}$ (F-49)

where zR has been omitted for simplicity. Equation (F-49) states that the gather Pp(h, ω) with water-bottom multiple attenuation is computed by taking the difference of the input gather P0(h, ω) and the input gather scaled by the water-bottom reflectivity rW(ω) and upward extrapolated by a depth interval twice the water depth zW.

Equation (F-49) applies to data for which CMP assumptions are acceptable. For complex water-bottom multiples, the horizontal layer assumption no longer holds and therefore input data must be the entire set of shot records along the seismic line [29]. As for free-surface multiple attenuation, however, computational efficiency in practice usually mandates use of individual CMP gathers in equation (F-49).

In practice, there are several limitations of the water-bottom multiple attenuation technique discussed here. First is the problem of an accurate estimate of the water-bottom reflectivity rW(ω) that is needed in equation (F-49). While attempts have been made to obtain an accurate and stable estimate of the water-bottom reflectivity [29], the problem remains as a major obstacle to the success of the technique in practice.

Additional issues arise in relation to spatial aliasing of the CMP data and missing near-offset data. As for the free-surface multiple attenuation, these problems may be circumvented by trace interpolation and extrapolation to zero offset using the Radon transform.

Perhaps the most serious limitation of the technique is the computational requirements to handle complex water-bottom geometry. In that case, the prediction and subtraction of multiples cannot be performed by the use of equation (F-49) based on phase-shift extrapolation. Instead, equation (F-45) needs to be employed wherein wave extrapolation has to be performed using the Kirchhoff integral to handle the irregular water-bottom geometry (layer replacement). The water-bottom geometry itself can be determined by migrating the stacked section using the constant-velocity Stolt migration (migration principles) and converting the picked time horizon associated with the water-bottom reflection to depth as for the layer replacement (layer replacement).

Figure F-12 shows a field data example of attenuation of complex water-bottom multiples [29]. The shot record in Figure F-12a exhibits strong multiples associated with a dipping water bottom. After the application of the technique based on equation (F-49), water-bottom multiples have been largely attenuated as shown in Figure F-12b. The stacked section shown in Figure F-12c associated with the shot record as in Figure F-12a contains complex water-bottom multiples. Any remaining multiples after multiple attenuation (Figure F-12d) may be attributed to the limitations in an accurate estimate of the water-bottom reflectivity and the adverse effect of the missing near-offset data.

### F.6 Spatial prediction filter

Consider a CMP-stacked data set P(x, t), where x is the CMP axis and t is the two-way zero-offset time axis. Apply Fourier transform in the t direction to decompose this 2-D data set to its frequency components P(x, ω). For each frequency component, define a complex vector P : P(x, ω) in the x direction.

We want to design a complex Wiener prediction filter F : F(x), such that, when applied to the input data vector P : P(x, ω), it yields an estimate of the input vector D : P(x + α, ω), at x + α, where D is the desired output vector and α is the prediction lag. For the specific problem of attenuating random noise uncorrelated from trace to trace, we shall consider the unit prediction lag α = 1. The output of the spatial prediction filtering is the predictable component of the input data vector that can be treated as the coherent signal, and the error in the prediction process is the random noise contained in the input data [30].

The prediction filtering is expressed by the following convolutional relation:

 ${\displaystyle Y(x)=P(x,\omega )\ast F(x),}$ (F-50)

where Y(x) represents the actual output from prediction filtering. We want to compute a complex prediction filter F(x) such that the difference E(x) between the actual output Y(x) and the desired output D(x) is minimum in the least-squares sense. Consider the discrete form of equation (F-50), with F(x) represented by the m–length complex vector F : (F0, F1, F2, …, Fm−1), P(x, ω) represented by the n–length complex vector P : (P0, P1, P2, …, Pn−1), and Y(x) represented by the (m + n − 1)-length complex vector Y : (Y0, Y1, Y2, …, Ym+n−1). Equation (F-50) can then be expressed in the following matrix form:

 ${\displaystyle {\begin{pmatrix}Y_{0}\\Y_{1}\\Y_{2}\\\vdots \\Y_{m+n-1}\end{pmatrix}}={\begin{pmatrix}P_{0}&0&\ldots &0\\P_{1}&P_{0}&\ldots &0\\P_{2}&P_{1}&P_{0}&\vdots \\\vdots &P_{2}&P_{1}&P_{0}\\P_{m-1}&\vdots &P_{2}&P_{1}\\0&P_{m-1}&\vdots &P_{2}\\\vdots &\ldots &P_{m-1}&\vdots \\0&0&\ldots &P_{m-1}\end{pmatrix}}{\begin{pmatrix}F_{0}\\F_{1}\\F_{2}\\\vdots \\F_{n-1}\end{pmatrix}}.}$ (F-51)

This is the equation for complete transient convolution. Define the (m + n − 1) × (m + n − 1) coefficient matrix on the right-hand side by L. Equation (F-51) takes the compact form

 ${\displaystyle \mathbf {Y} =\mathbf {LF} .}$ (F-52)

The error vector E : (E0, E1, E2, …, Em+n−1) is defined as the difference between the desired output D : (D0, D1, D2, …, Dm+n−1) and the actual output Y : (Y0, Y1, Y2, …, Ym+n−1):

 ${\displaystyle \mathbf {E} =\mathbf {D} -\mathbf {Y} .}$ (F-53)

By substituting equation (F-52) into equation (F-53), we obtain

 ${\displaystyle \mathbf {E} =\mathbf {D} -\mathbf {LF} .}$ (F-54)
Figure F-12  (a) A shot record with strong water-bottom multiples, (b) the same record after water-bottom multiple attenuation, (c) a portion of the stacked section associated with the data as in (a) that exhibits complex water-bottom multiples, (d) a portion of the stacked section associated with the data as in (b) after water-bottom multiple attenuation. [29].

The energy of the error vector is

 ${\displaystyle S=\mathbf {E^{T\ast }E} ,}$ (F-55)

where T denotes matrix transpose and * denotes complex conjugate. Now, substitution of equation (F-54) into the definition defined by equation (F-55) yields:

 ${\displaystyle S=(\mathbf {D} -\mathbf {LF} )^{\mathbf {T\ast } }(\mathbf {D-LF} )}$ (F-56)

Expand the right-hand side:

 ${\displaystyle S=\mathbf {D^{T\ast }D} -\mathbf {F^{T\ast }L^{T\ast }D-D^{T\ast }LF+F^{T\ast }L^{T\ast }LF} .}$ (F-57)

We want to estimate a prediction filter vector F such that the quantity S is minimum. This condition leads to setting the derivative of S with respect to F to zero. Differentiate both sides of equation (F-57) with respect to F and observe the requirement for least-squares minimization that ∂S/F = 0:

 ${\displaystyle -\mathbf {D^{T\ast }L+F^{T\ast }L^{T\ast }L=0} .}$ (F-58)

Apply matrix transpose and rearrange the terms:

 ${\displaystyle (\mathbf {L^{T\ast }L)^{T\ast }F} =\mathbf {L^{T\ast }D} .}$ (F-59)

Now define

 ${\displaystyle \mathbf {G} =\mathbf {L^{T\ast }D} }$ (F-60a)

and

 ${\displaystyle \mathbf {R} =\mathbf {L^{T\ast }L} ,}$ (F-60b)

which yields the relation

 ${\displaystyle \mathbf {R^{T\ast }} =\mathbf {R} ,}$ (F-60c)

thus making R a Hermitian matrix of the size (m + n − 1) × (m + n − 1).

Use the relations given by equations (F-60a,F-60b,F-60c), to rewrite equation (F-58) in the form

 ${\displaystyle \mathbf {RF} =\mathbf {G} ,}$ (F-61)

which can be rewritten by using equations (F-60a, F-60b) as follows:

 ${\displaystyle \mathbf {F} =\mathbf {(L^{T\ast }L)^{-1}L^{T\ast }D} .}$ (F-62)

This solution is of the same form as the generalized linear inverse form of the discrete Radon transform represented by equation (F-20). As for the latter application, singular-value decomposition technique based on equation (F-23) can be used to solve for the complex Wiener filter coefficients F : (F0, F1, F2, …, Fm−1).

An efficient recursive algorithm to solve for the complex Wiener prediction filter coefficients F in equation (F-61) is described by Treitel [31]. The technique makes use of the Hermitian property (equation F-60c) of the autocorrelation matrix R. Write this matrix in terms of its real Rr and imaginary Ri components as

 ${\displaystyle \mathbf {R} =\mathbf {R} _{r}+i\mathbf {R} _{i}.}$ (F-63a)

Similarly, write the complex Wiener filter F and the crosscorrelation matrix G in the same form:

 ${\displaystyle \mathbf {F} =\mathbf {F} _{r}+i\mathbf {F} _{i}}$ (F-63b)

and

 ${\displaystyle \mathbf {G} =\mathbf {G} _{r}+i\mathbf {G} _{i}.}$ (F-63c)

Substitute equations (F-63a,F-63b,F-63c) into equation (F-61) and equate the real and imaginary parts of boths sides of the resulting expression to get

 ${\displaystyle \mathbf {R} _{r}\mathbf {F} _{r}-\mathbf {R} _{i}\mathbf {F} _{i}=\mathbf {G} _{r}}$ (F-64a)

and

 ${\displaystyle \mathbf {R} _{i}\mathbf {F} _{r}+\mathbf {R} _{r}\mathbf {F} _{i}=\mathbf {G} _{i}.}$ (F-64b)

Finally, write equations (F-64a, F-64b) in matrix form:

 ${\displaystyle {\begin{pmatrix}\mathbf {R} _{r}&-\mathbf {R} _{i}\\\mathbf {R} _{i}&\mathbf {R} _{r}\end{pmatrix}}{\begin{pmatrix}\mathbf {F} _{r}\\\mathbf {F} _{i}\end{pmatrix}}={\begin{pmatrix}\mathbf {G} _{r}\\\mathbf {G} _{i}\end{pmatrix}}.}$ (F-65)

The square matrix on the left-hand side of this equation is a block Toeplitz matrix [32] [33]. This property leads to an efficient recursive algorithm that does not require complex arithmetic to solve for the complex prediction filter coefficients F [31].

## References

1. Pekeris, 1948, Pekeris, C. L., 1948, Theory of propagation of explosive sounds in shallow water: Geol. Soc. Am. Mem. 27.
2. Press and Ewing, 1950, Press, F. and Ewing, M., 1950, Propagation of explosive sound in a liquid layer overlying a semi-infinite elastic solid: Geophysics, 15, 111–148.
3. Ewing et al. (1957), Ewing, M., Jardetsky, W. S. and Press, F., 1957, Elastic waves in layered media: McGraw-Hill Book Co.
4. McMechan and Yedlin (1981), McMechan, G. and Yedlin, M., 1981, Analysis of dispersive waves by wavefield transformation: Geophysics, 46, 869–879.
5. Yilmaz (1981), Yilmaz, O., 1981, Modeling guided waves in an irregular water layer: Presented at the 51st Ann. Internat. Mtg. Soc. Expl. Geophys.
6. Ottolini, 1982, Ottolini, R., 1982, Migration of seismic data in angle-midpoint coordinates: Ph.D. thesis, Stanford University.
7. Clayton and McMechan (1981), Clayton, R. W. and McMechan, G., 1981, Inversion of refraction data by wave-field extrapolation: Geophysics, 46, 860–868.
8. Beylkin, 1987, Beylkin, G., 1987, Discrete Radon transform: IEEE Trans. on Acoustics, Speech and Signal Processing, AASP-35, 162–172.
9. Thorson and Claerbout (1985), Thorson, J. R. and Claerbout, J. F., 1985, Velocity-stack and slant-stack stochastic inversion: Geophysics, 50, 2727–2741.
10. Yilmaz, 1989, Yilmaz, O., 1989, Velocity-stack processing: Geophys. Prosp., 37, 357–382.
11. Foster and Mosher, 1992, Foster, D. J. and Mosher, C. C., 1992, Suppression of multiple reflections using the Radon transform: geophysics, 57, 386–395.
12. Hampson (1986), Hampson, D., 1986, Inverse velocity stacking for multiple elimination: J. Can. Soc. Expl. Geophys., 22, 44–55.
13. 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.
14. Press et al., 1986, Press, W. H., Flannery, B. P., Teukolsky, S. A. and Vetterling, W. T., 1986, Numerical Recipes: The Art of Scientific Computing: Cambridge University Press.
15. Kostov, 1990, Kostov, C., 1990, Toeplitz structure in slant-stack inversion: 60th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 1618–1621.
16. Hampson, 1987, Hampson, D., 1987, The discrete Radon transform: A new tool for image enhancement and noise suppression: 57th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 141–143.
17. Riley and Claerbout, 1976, Riley, D. C. and J. F. Claerbout, 1976, 2-D multiple reflections: Geophysics, 41, 592–620.
18. Kennett (1979), Kennett, B. L. N., 1979, The suppression of surface multiples on seismic records: Geophys. Prosp., 28, 584–600.
19. Verschuur (1991), Verschuur, D. J., 1991, Surface-related multiple elimination, an inversion approach: Ph. D. thesis, Delft Univ. of Tech.
20. Verschuur et al. (1992), Verschuur, D. J., Berkhout, A. J. and Wapenaar, C. P. A., 1992, Adaptive surface-related multiple elimination: Geophysics, 57, 1166–1177.
21. Carvalho et al. (1991), Carvalho, P. M., Weglein, A. B. and Stolt, R. H., 1991, Examples of a nonlinear inversion method based on the T matrix of scattering theory: Application to multiple suppression: 61st Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 1319–1322.
22. Dragoset and McKay (1993), Dragoset, W. H. and MacKay, S., 1993, Surface multiple elimination and subsalt imaging: Expl. Geophys., 24, 463–472.
23. Verschuur and Berkhout (1994), Verschuur, D. J. and Berkhout, A. J., 1994, Multiple technology, part 1: Estimation of multiple reflections: 64th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 1493–1496.
24. Berkhout and Verschuur (1995), Berkhout, A. J. and Verschuur, E. J., 1995, Estimation of multiple scattering by iterative inversion, Part 1: Theoretical considerations: 55th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 715–718.
25. Verschuur et al. (1995), Verschuur, D. J., Berkhout A. J., and Kelamis, P. G., 1995, Estimation of multiple scattering by iterative inversion — part II: examples of marine and land data: 65th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 1470–1473.
26. Kelamis and Verschuur (1996), Kelamis, P. G. and Verschuur, D. J., 1996, Multiple elimination strategies for land data: Eur. Assoc. Expl. Geophys. Mtg., Extended Abstracts.
27. Dragoset and Jericevic (1998), Dragoset, W. H. and Jericevic, Z., 1998, Some remarks on surface multiple elimination: Geophysics, 63, 772–789.
28. Morley (1982), Morley, L., 1982, Predictive techniques for marine multiple suppression: Ph. D. dissertation, Stanford Univ.
29. Wiggins (1988), Wiggins, J. W., 1988, Attenuation of complex water-bottom multiple by wave-equation-based prediction and subtraction: Geophysics, 53, 1527–1539.
30. Canales (1984), Canales, L., 1984, Random noise reduction: 54th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 525.
31. Treitel (1974), Treitel, S., 1974, The complex Wiener filter: Geophysics, 39, 169–173.
32. Robinson, 1967, Robinson, E. A., 1967, Multichannel time series analysis with digital computer programs: Holden-Day.
33. Robinson and Treitel, 1980, Robinson, E. A. and Treitel, S., 1980, Geophysical signal analysis: Prentice-Hall, Inc.