# Topics in moveout and statics corrections

Series | Investigations in Geophysics |
---|---|

Author | Öz Yilmaz |

DOI | http://dx.doi.org/10.1190/1.9781560801580 |

ISBN | ISBN 978-1-56080-094-1 |

Store | SEG Online Store |

## Contents

- 1 C.1 The shifted hyperbola
- 2 C.2 Moveout stretch
- 3 C.3 Equations for a dipping reflector
- 4 C.4 Traveltime decomposition for residual statics estimation
- 5 C.5 Depth estimation from refracted arrivals
- 6 C.6 Equations for a dipping refractor
- 7 C.7 The plus-minus times
- 8 C.8 Generalized linear inversion of refracted arrivals
- 9 C.9 Refraction traveltime tomography
- 10 C.10
*L*_{1}-norm refraction statics - 11 Figures
- 12 Equations
- 13 See also
- 14 References
- 15 External links

## C.1 The shifted hyperbola

The objective in this section is to review the higher-order accuracy in normal moveout for a horizontally layered earth model. Refer to the traveltime equation (**3**) by Tanner and Koehler, 1969^{[1]}:

**(**)

where *x* is the offset, and

**(**)

where *t*_{0} is the two-way zero-offset traveltime

**(**)

**(**)

and

**(**)

with the additional definitions of the terms ^{[2]}

**(**)

Note that

**(**)

since

**(**)

where Δ*t _{i}* is the vertical two-way traveltime through the

*i*-th layer,

*v*is the velocity, and Derivation of equation (

**C-1**) is based on the parametric equations for traveltime and offset for a horizontally layered earth model

^{[3]};

^{[4]}and is given by

^{[2]}.

Drop higher-order terms in equation (**C-1**) to obtain the *fourth-order moveout* equation

**(**)

Now, substitute for the coefficients from equations (**C-2a**), (**C-2b**) and (**C-3a**) to obtain

**(**)

This equation can in principle be used to compute a velocity spectrum (Velocity analysis). First, drop the fourth-order term to get the small-spread hyperbolic equation

**(**)

Compute the velocity spectrum using equation (**C-6**), and pick an initial velocity function *v _{rms}*(

*t*

_{0}). Then, use this picked velocity function in equation (

**C-5**) and compute a velocity spectrum for the parameter

*C*

_{2}. Pick a function

*C*

_{2}(

*t*

_{0}) and use in equation (

**C-5**) to recompute the velocity spectrum. Finally, pick an updated velocity function

*v*(

_{rms}*t*

_{0}) from this velocity spectrum.

The fourth-order moveout equation (**C-4**) can be expressed exactly by a time-shifted hyperbolic traveltime equation of the form ^{[2]}

**(**)

where *S* is a constant of the form

**(**)

For *S* = 1, equation (**C-7**) reduces to the conventional moveout equation (**C-6**).

Figure C-1 shows the traveltime trajectories of the hyperbolic moveout equation (**C-6**) and the time-shifted hyperbolic equation (**C-7**) where *τ _{s}* =

*t*

_{0}(1 − 1/

*S*). Note that the shifted hyperbola is a better match to the true traveltime trajectory at far offsets. The latter was computed for a horizontally layered earth model using the exact parametric equations for traveltime and offset

^{[3]};

^{[4]}.

The shifted hyperbola equation satisfies the requirements for a moveout equation:

- The traveltime trajectory is symmetric with respect to the time axis.
- The instantaneous velocity
*dx*/*dt*is never zero. - For constant velocity, it reduces to the exact hyperbolic moveout equation.

As for the fourth-order moveout equation (**C-5**), this equation can, in principle, be used to conduct the velocity analysis of CMP gathers. First, set *S* = 1 in equation (**C-7**) to get equation (**C-6**). Compute the velocity spectrum using equation (**C-6**), and pick an initial velocity function *v _{rms}*(

*t*

_{0}). Then, use this picked velocity function in equation (

**C-7**) and compute a velocity spectrum for the parameter

*S*. Pick a function

*S*(

*t*

_{0}) and use in equation (

**C-7**) to recompute the velocity spectrum. Finally, pick an updated velocity function

*v*(

_{rms}*t*

_{0}) from this velocity spectrum.

De Bazelaire^{[5]} offers an alternative moveout equation to achieve higher-order accuracy at far offsets:

**(**)

where *t*_{0} is the two-way zero-offset time, *t _{p}* is related to the time at which the asymptotes of the hyperbolic traveltime trajectory converge (Figure C-2), and

*v*is the reference velocity assigned to the layer below the recording surface (not the near-surface layer). When

_{s}*t*=

_{p}*t*

_{0}, equation (

**C-9**) reduces to the small-spread hyperbolic equation (

**C-6**).

## C.2 Moveout stretch

Refer to Figure 3.1-10 for a sketch of a wavelet before and after moveout correction. The moveout equation is

**(**)

where *t* is the two-way traveltime associated with a source-receiver separation *x*, *t*_{0} is the two-way vertical traveltime — the time after moveout correction, and v is the moveout velocity.

Consider a reflection event represented by a wavelet of dominant period *T* with an arrival time *t* at offset *x*. After normal-moveout correction, the dominant period becomes *T*_{0} = *T* + Δ*T*. The moveout equation (**C-10**) is associated with the onset of the wavelet. Similarly, the moveout equation with the termination of the wavelet is expressed by

**(**)

Expand the terms on both sides:

**(**)

By making the substitution from equation (**C-10**), we obtain

**(**)

Simplify and rearrange the terms

**(**)

Now, ignore the second term on the right-hand side of the equation and observe that Δ*t _{NMO}* =

*t*−

*t*

_{0}to obtain

**(**)

Assume that *t*_{0} >> *T* and rearrange the terms to obtain a relationship for change in the period of the wavelet as a result of moveout correction:

**(**)

Now, we want to express equation (**C-13**) in terms of the dominant frequency *f* of the wavelet. Start with the relation

**(**)

and obtain

**(**)

Finally, combine equations (**C-13**) and (**C-15**) to obtain the equation for the absolute value of frequency stretching:

**(**)

This is the same as equation (**6**) in the main text.

## C.3 Equations for a dipping reflector

We want to derive equation (**7**) of the main text using the geometry of a dipping reflector as shown in Figure C-3. (The derivation provided here is courtesy Zhiming Li, 1999). The distance from source location *S* to reflection point *R* back to receiver location *G* is given by *SR* + *RG* = *FR* + *RG* = *FG* = *vt*, where point *F* is the mirror image of the source location *S*, *v* is the velocity of the medium above the dipping reflector, and *t* is the traveltime from *F* to *G*.

To compute the traveltime *t* associated with the distance *FG*, we shall need to compute the coordinates of *F* : (*x _{F}*,

*z*) and

_{F}*G*: (

*x*,

_{G}*z*), so that

_{G}

**(**)

From the geometry of Figure C-3, the coordinates of *F* : (*x _{F}*,

*z*) are

_{F}and

where *ϕ* is the dip angle of the reflector. Substitute the relation *OF* = *OS* = *OM* + *x*/2 into the above equations to obtain

**(**)

and

**(**)

Again, from the geometry of Figure C-3, the coordinates of *G* : (*x _{G}*,

*z*) are

_{G}

**(**)

and

**(**)

Substitute equations (**C-18a**,**C-19a**,**C-17**)

**(**)

then use the trigonometric relation cos 2*ϕ* = 2 cos^{2} *ϕ* − 1 and apply some algebra to obtain

**(**)

Further algebraic manipulation yields

**(**)

Finally, note from the geometry of Figure C-3 the relation

**(**)

Subsitute equation (**C-22**) into equation (**C-21**), and use the relations 2*MN* = *vt*_{0} and *FG* = *vt* to obtain the desired expression for the reflection traveltime associated with a dipping reflector

**(**)

and the moveout velocity

**(**)

Equations (**10**) and (**11**) of the main text.

## C.4 Traveltime decomposition for residual statics estimation

We want to model traveltime deviations associated with a reflection event on moveout-corrected CMP gathers by the following equation ^{[6]}; ^{[7]}:

**(**)

where *s _{j}* is the residual statics shift at the

*j*th source location,

*r*is the residual statics shift at the

_{i}*i*th receiver location,

*G*is the structure term at the

_{k}*k*th midpoint location, [

*k*= (

*i*+

*j*)/2], is the residual moveout at the

*k*th midpoint location.

For *m* picks of *t _{ij}*, and

*n*shot locations,

_{s}*n*receiver locations,

_{r}*n*midpoint locations, we have the following set of equations:

_{G}

**(**)

A more parsimonious parameterization can be achieved by an appropriate traveltime picking strategy such that the structure term can be set to zero. Additionally, the residual moveout term coefficient *M* may be set to a constant. The model equation (**C-24a**) can be rewritten as

**(**)

Equation (**C-24b**) is modified, accordingly

**(**)

Whichever the preferred modeled traveltimes, write the respective equations (**C-24b**) or (**C-25b**) in matrix notation as

**(**)

where **t′** is the column vector of *m*-length (number of traveltime picks) in equation (**C-24b**) or (**C-25b**), **L** is the sparse matrix in the same equations with dimensions *m* × (*n _{s}* +

*n*+

_{r}*n*+

_{G}*n*) in case of equation (

_{G}**C-24b**) and

*m*× (

*n*+

_{s}*n*+ 1) in case of equation (

_{r}**C-25b**), and

**p**is the column vector of (

*n*+

_{s}*n*+

_{r}*n*+

_{G}*n*)-length in case of equation (

_{G}**C-24b**) and (

*n*+

_{s}*n*+ 1)-length in case of equation (

_{r}**C-25b**) on the right-hand side of these equations. Except the three elements in each row, the

**L**matrix contains zeros.

The generalized least-squares solution to equation (**C-26**) satisfies the requirement that the energy of the error vector

**(**)

is minimum ^{[8]}. The cumulative error energy is

**(**)

Substitute equation (**C-27a**) into equation (**C-27b**) and use equation (**C-26**) to obtain

**(**)

Minimization of *E* with respect to *p* requires that

**(**)

which, in the case of equation (**C-24b**), yields *n _{s}* +

*n*+

_{r}*n*+

_{G}*n*equations and that many unknowns. These equations can be solved for the residual statics associated with

_{G}*n*source locations,

_{s}*n*receiver locations,

_{r}*n*structural terms, and

_{G}*n*residual moveout terms.

_{G}Finally, the solution that satisfies the minimization requirement given by equation (**C-27d**) follows ^{[8]}

**(**)

where **t** denotes the column vector of *m*-length that represents the traveltime deviations picked from moveout-corrected CMP gathers, and *T* denotes the matrix transposition.

A practical scheme for solving equation (**C-24a**) is based on the Gauss-Seidel method. In this scheme, each term on the right-hand side of equation (**C-24a**) is computed by the set of recursive equations

**(**)

**(**)

**(**)

and

**(**)

where *i* and *j* are the receiver and source indexes, respectively, *l* = |*i − j*| is the offset index, *k* = (*i* + *j*)/2 is the midpoint index, *m* is the iteration index, and *n _{h}* is the fold of coverage. The solutions in equations (

**C-29**) are based on the orthogonality of the shot and receiver axes, and the orthogonality of the midpoint and offset axes. Equations (

**C-29**) can be modified as

**(**)

**(**)

**(**)

and

**(**)

This modification enables us to compute and store the sum of the traveltime deviations ∑*t _{ij}* thus circumventing the need for storing the individual traveltime deviations. The process is iterated until an index

*m*that yields the least-squares solution such that the differences between the solutions for

*s*,

_{j}*r*,

_{i}*G*, and

_{k}*M*from the

_{k}*m*th and (

*m*− 1)st iteration are smaller than some specified values.

## C.5 Depth estimation from refracted arrivals

Refer to Figure 3.4-11a for a sketch of refracted arrivals associated with a flat refractor. We want to estimate the depth to bedrock *z _{w}* (thickness of the weathering layer) at some shot-receiver station. The traveltime

*t*for the critically refracted arrival associated with a shot

*S*and receiver

*R*has three parts

**(**)

The terms in this equation can be expressed in terms of the near-surface model parameters

**(**)

where *x* is the shot-receiver separation, *z _{w}* is the depth to bedrock, and

*v*and

_{w}*v*are weathering and bedrock velocities, respectively. The critical angle of refraction

_{b}*θ*is related to the velocities by the relation

_{c}

**(**)

Combine the first and third terms on the right, and split the second term in equation (**C-31b**) to get

**(**)

Substitute equation (**C-32**) into equation (**C-33a**) and apply some algebra to obtain

**(**)

Note that this is the equation of a line

**(**)

with its slope given by 1/*v _{b}* — inverse of the bedrock velocity, and the intercept time

*t*given by

_{i}

**(**)

By measuring the slope of the line associated with refracted arrivals in Figure 3.4-11a, we estimate the bedrock velocity *v _{b}*. We also estimate the intercept time

*t*by extending the line to

_{i}*x*= 0. Finally, we assume a value for the weathering velocity

*v*. Then, equation (

_{w}**C-35**) can be rearranged to compute the depth to bedrock

*z*as

_{w}

**(**)

This is equation (**41a**) in the text.

The depth to bedrock can also be computed from the crossover distance *x _{c}* where the refracted arrival and the direct arrival coincide (Figure 3.4-11a). The equation for the direct arrivals is

**(**)

whereas equation (**C-34**) describes the refracted arrivals. At the crossover distance *x _{c}*, we equate the arrival times in equations (

**C-34**) and (

**C-37**) as

**(**)

By substituting equation (**C-35**), we obtain

**(**)

By simplifying, we obtain the expression for depth to bedrock in terms of the crossover distance:

**(**)

This is equation (**41b**) in the text.

## C.6 Equations for a dipping refractor

Consider the dipping refractor depicted in Figure 3.4-11c. The traveltime for the refracted arrivals from downdip shooting is given by

**(**)

From the geometry and the variables of Figure 3.4-11c, the terms in equation (**C-40a**) are explicitly written as

**(**)

Apply some algebra to obtain

**(**)

This is the equation of a line

**(**)

with its inverse slope defined as

**(**)

and the intercept time given by

**(**)

Equation (**C-41a**) can be rewritten for the refracted arrivals from updip shooting (Figure 3.4-11d) as

**(**)

Again, this is the equation of a line

**(**)

with its inverse slope defined as

**(**)

and the intercept time given by

**(**)

To compute the refractor dip and bedrock velocity, first, rewrite equations (**C-41c**) and (**C-42c**)

**(**)

and

**(**)

By subtracting equation (**C-43b**) from (**C-43a**), we obtain an expression for the refractor dip *φ*:

**(**)

This is equation (**46a**) of the text. Measure the slopes of the downdip and updip refracted arrivals in Figure 3.4-11d — respectively. Then, assume a value for the weathering velocity *v _{w}*. By direct substitution into equation (

**C-44**), compute the refractor dip

*φ*.

To obtain the bedrock velocity *v _{b}*, first, rewrite slopes from equations (

**C-41c**) and (

**C-42c**) as

**(**)

and

**(**)

Add equations (**C-45a**) and (**C-45b**) and take the inverse. Then use the relation given by equation (**C-32**) to obtain

**(**)

This is equation (**46b**) of the text. Finally, we compute the depth to the bedrock at shot-receiver stations by substituting the estimates for refractor dip *φ* from equation (**C-44**) and the bedrock velocity *v _{b}* from equation (

**C-46**), and the measurements for the intercept times from the arrival times in Figure 3.4-11d into equations (

**C-41d**) and (

**C-42d**) to obtain

**(**)

and

**(**)

These equations are equivalent to equation (**46c**) of the text.

## C.7 The plus-minus times

Consider the three raypaths in Figure 3.4-12a associated with shot-receiver pairs *AD*, *DG*, and *AG*. The plus and minus times are defined as

**(**)

and

**(**)

The times given on the right side of these equations are the picked values from the first breaks for the three raypaths shown in Figure 3.4-12a. From the raypath configuration, we have the relation

**(**)

By using the geometric relations from Figure (3.4-12a), equation (**C-49a**) can be expressed in terms of the critical angle of refraction *θ* and depth to bedrock *z _{w}* as

**(**)

Finally, by using the relation (**C-32**), we have the relation for the plus time *t _{+}* in terms of the near-surface parameters — depth to bedrock

*z*, weathering velocity

_{w}*v*, and bedrock velocity

_{w}*v*:

_{b}

**(**)

This is equation (**48c**) of the main text.

Now, consider the minus time *t _{_}* defined by equation (

**C-48b**). By using the raypath configuration depicted in Figure 3.4-12b, we have

**(**)

Make substitutions in terms of the near-surface parameters to get

**(**)

By some algebraic manipulation, we get

**(**)

where *x* is the source-receiver separation *AD*. Compare the first two terms on the right with equation (**C-48b**), and rewrite equation (**C-49c**) to obtain

**(**)

This is equation (**48a**) of the text.

## C.8 Generalized linear inversion of refracted arrivals

For an arbitrary source-receiver geometry, given a set of observed traveltimes associated with the refracted arrivals, we can estimate the parameters associated with a single-layer near-surface model using the generalized linear inversion (GLI). The parameters consist of the refractor velocity, and the velocity and thickness of the near-surface layer at all shot/receiver locations. The GLI solution for the parameters satisfies the requirement that the difference between the observed (picked) refracted arrival times and the estimated (modeled) times is minimum in the least-squares sense. The modeled times are computed using the traveltime equation for refracted arrivals for a flat refractor considered as the base of a weathering layer.

The GLI schemes that allow velocity and thickness of a near-surface layer to vary spatially require iterative strategies ^{[9]}^{[10]}^{[11]}. Starting with initial estimates for the near-surface layer velocity and thickness, and an initial estimate for the bedrock velocity, these parameters are changed such that the difference between the observed (picked) refracted arrival times and the estimated (modeled) times is minimum in the least-squares sense. The GLI method is not only applicable to 2-D line recording but also to 3-D swath recording geometries ^{[12]}^{[13]}. It is important to parameterize the near-surface layer parsimoniously while conforming with the basic assumptions required for the use of refracted arrivals in estimating the near-surface model.

There are several ways to parameterize the near-surface layer. The most general formulation would include varying the weathering and the bedrock velocities and the thickness of the weathering layer at all shot/receiver locations. This, however, would require linearizing the problem and iterating over the estimated parameters ^{[9]}^{[10]}^{[11]}. In a simplified version of this general formulation, the weathering velocity may be fixed and is assumed to be known. This leaves the weathering thickness and bedrock velocity as spatially varying parameters. A further simplification may be made by defining a polynomial as a function of the space variable for the bedrock velocity ^{[14]}. In practice, we often find that iterative GLI schemes suffer from stability problems. When an inversion scheme leads to unstable solutions, the most likely reason is the ill-posed nature of the problem.

In this section, two robust parameterizations of the near-surface layer are presented and the GLI method is used to estimate the associated parameters. In each, we assume that the near-surface is made up of a single weathering layer with a significant velocity contrast at its base. The *variable-thickness* scheme allows the thickness of the weathering layer to vary spatially, while assuming a fixed known value for the weathering velocity. The bedrock velocity, however, is included in the parameterization and is assumed to be constant. The *variable-velocity* scheme allows the weathering velocity to vary spatially, while fixing the refractor position at a specified depth. The bedrock velocity, again, is treated as a parameter to be estimated and is asssumed to be constant. We find that in many field data applications, these two inversion schemes are able to remove long-wavelength statics variations from the data.

We want to describe the near-surface with minimal parameterization and consider the model shown in Figure 3.4-13. The traveltime for the refracted raypath from the shot location *S _{j}* to the receiver location

*R*is given by

_{i}

**(**)

The first and the third terms are associated with the raypaths within the weathering layer and the second term is associated with the raypath within the bedrock along the refractor. In Figure 3.4-13, *θ _{c}* is the critical angle of refraction which is expressed in terms of the weathering and bedrock velocities by the relation

*θ*= sin

_{c}^{−1}(

*v*/

_{w}*v*). Also, as depicted in Figure 3.4-13, we assume a flat refractor. When refractor dip is taken into consideration, the problem cannot be readily linearized.

_{b}By rewriting equation (**C-50**), for a flat or near-flat refractor, we obtain

**(**)

By regrouping the terms, we get

**(**)

Finally, by rewriting equation (**C-52**) in terms of the near-surface parameters, we obtain the model equation for the refracted arrivals:

**(**)

In additon to asumming a flat refractor, we fix the bedrock velocity but retain it as a parameter to be estimated. Under these assumptions, equation (**C-52**) can be rewritten in the following form:

**(**)

where

**(**)

**(**)

and

**(**)

Note that *T _{j}* and

*T*are the intercept time anomalies at shot and receiver locations, respectively, and

_{i}*s*is the bedrock slowness. Hence, for

_{b}*n*shot/receiver stations the parameter vector is

**p**: (

*T*

_{1},

*T*

_{2}, …,

*T*;

_{n}*s*). We refer to the scheme based on equation (

_{b}**C-54**) as the

*variable-thickness*GLI solution. Once the parameter vector

**p**: (

*T*

_{1},

*T*

_{2}, …,

*T*;

_{n}*s*) is estimated, then the thickness of the weathering layer below shot and receiver locations can be computed using equations (

_{b}**C-55**) and (

**C-56**), respectively. For

*m*picks of

*t*, and

_{ij}*n*+ 1 parameters

**p**: (

*T*

_{1},

*T*

_{2}, …,

*T*;

_{n}*s*), we have the following set of equations:

_{b}

**(**)

By writing these equations in matrix notation, we have

**(**)

where **t′** is the column vector of *m*-length on the left-hand side of equation (**C-58**), **L** is the sparse matrix in the same equation with dimensions *m* × (*n* +1), and **p** is the column vector of (*n* + 1)-length on the right-hand side of equation (**C-58**). Except for the three elements in each row, the **L** matrix contains zeros. The GLI solution to equation (**C-59**) satisfies the requirement that the energy of the error vector

**(**)

is minimum and is given by:

**(**)

where **t** denotes the column vector of *m*-length that represents the observed (picked) refracted arrival times, and *T* denotes matrix transposition.

In summary, the variable-thickness scheme for refraction statics is as follows:

- Assume a value for the weathering velocity
*v*. This can be varied spatially based on available uphole information._{w} - Estimate the parameter vector
**p**: (*T*_{1},*T*_{2}, …,*T*;_{n}*s*), hence compute the intercept time anomalies at shot/receiver locations and the bedrock slowness by solving equation (_{b}**C-61**). - Solve equations (
**C-55**) and (**C-56**) for the weathering layer thickness at shot and receiver stations, respectively.

The shape of the weathering layer derived from the variable-thickness scheme strictly depends on the assumed value for the weathering velocity. To demonstrate this important aspect of the variable-thickness scheme, consider the near-surface model in Figure C-4a with a flat refractor R1, and constant weathering and bedrock velocities. We use equation (**C-53**) and compute the refracted arrival times associated with shot/receiver locations as indicated in Figure C-4a. Then we use these arrival times in equation (**C-61**) and assume a value for the weathering velocity different from its true value to estimate the thickness of the weathering layer at all shot/receiver stations. Figures C-4b,c show the estimated refractor shapes R2, R3 for two different weathering velocities. Note the significant departure from the true refractor shape R1. Results of Figure C-4 clearly demonstrate that the estimated refractor shape using the variable-thickness scheme (whether it is based on GRM or GLI) does not yield the true refractor shape. Instead, the uncertainty in weathering velocity significantly influences the implied refractor shape.

We now consider an alternative parameterization for the near-surface model. We assume a flat base of weathering as in Figure 3.4-13. This assumption makes the weathering thickness a known quantity in equation (**C-53**), and leaves the weathering and bedrock velocities as unknown parameters to be estimated. We take the similar view as in the variable-thickness scheme, and further assume the bedrock velocity to be a constant parameter. Under these implicit constraints, equation (**C-53**) takes the form

**(**)

where

**(**)

**(**)

**(**)

and

**(**)

**Figure C-4**(a) A single-layered near-surface model with a flat refractor R1;*v*= 3000 f/s,_{w}*v*= 9000 f/s. Shot locations are denoted by × and receiver locations are denoted by the vertical bars. Traveltimes computed from the model in (a) and equation (_{b}**C-53**) are put into the variable-thickness GLI inversion scheme (equation 2-5) to obtain the results shown in (b) where the weathering velocity is assumed to be 3500 f/s, and (c) where the weathering velocity is assumed to be 2500 f/s. The estimated refractor shapes are denoted by R2 and R3. Note the departure from the true location of the flat refractor R1 due to the uncertainty in the weathering velocity.**Figure C-5**(a) Same model as in Figure C-4a. (b) Results of the variable-velocity GLI inversion scheme (equation 2-13) using the traveltimes computed from equation (**C-53**) and the model in (a) with an assumed refractor position R2 that is different from the true position R1 of the refractor. The true weathering and bedrock velocities,*v*and_{w}*v*, are 3000 f/s and 9000 f/s, respectively. Shot locations are denoted by × and receiver locations are denoted by vertical bars. While the GLI inversion yields the true refractor velocity, note the spatially varying adjustment of the weathering velocity to the change of the refractor position from R1 to R2 to compensate for the thickness R1-R2._{b}

Hence, for *n* shot/receiver stations, the parameter vector is **p** : (*α*_{1}, *α*_{2}, …, *α _{n}*;

*s*). We refer to the scheme based on equation (

_{b}**C-62**) as the

*variable-velocity*GLI solution. Once the parameter vector

**p**: (

*α*

_{1},

*α*

_{2}, …,

*α*;

_{n}*s*) is estimated, then the weathering velocity below shot and receiver locations can be computed using equations (

_{b}**C-63**) and (

**C-64**), respectively. For

*m*picks of

*t*, and

_{ij}*n*+ 1 parameters

**p**: (

*α*

_{1},

*α*

_{2}, …,

*α*;

_{n}*s*), we have the following set of equations:

_{b}

**(**)

We write these equations in matrix notation as in equation (**C-59**), where **t′** is the column vector of *m*-length on the left-hand side of equation (**C-67**), **L** is the sparse matrix in the same equation with dimensions *m* × (*n* + 1), and **p** is the column vector of (*n* + 1)-length on the right-hand side of equation (**C-67**). Except for the three elements in each row, the **L** matrix contains zeros. The GLI solution is given by equation (**C-61**) where the **L** matrix is defined as in equation (**C-67**).

In summary, the variable-velocity scheme for refraction statics is as follows:

- Specify a flat datum to which shot and receivers are to be lowered along vertical raypaths.
- Estimate the parameters:
**p**: (*α*_{1},*α*_{2}, …,*α*;_{n}*s*) using the GLI solution given by equation (_{b}**C-61**) where the**L**matrix is defined as in equation (**C-67**). - Solve equations (
**C-63**) and (**C-64**) for the weathering velocity at shot and receiver stations, respectively.

The estimated weathering velocity from the variable-velocity scheme depends on the assumed refractor position. Consider the near-surface model in Figure C-5a, with a flat refractor R1, and constant weathering and bedrock velocities. We use equation (**C-53**) and compute the refracted arrival times associated with shot/receiver locations as indicated in Figure C-5a. We then use these arrival times in equation (**C-61**) and assume a refractor position R2 to estimate the weathering velocity at all shot/receiver stations. Figure C-5b shows the estimated weathering velocity which departs from the true value *v _{w}*. The difference arises from compensating for the difference between the true refractor position R1 and the assumed refractor position R2 (Figure C-5a).

In equations (**C-58**) and (**C-67**), we made no distinction between a shot and a receiver if they occupy the same location on the surface. This strictly surface-consistent solution is of course not valid for a dynamite source. However, shots can be brought up to surface with an uphole correction prior to setting up equations (**C-58**) and (**C-67**).

For most field data applications, both the variable-thickness and the variable-velocity schemes yield comparable statics solutions. The solution from the variable-thickness scheme is sensitive to the uncertainty in weathering velocity, and the solution from the variable-velocity scheme is influenced by the assumed depth of the flat refractor.

## C.9 Refraction traveltime tomography

Return to equation (**C-52**) and re-express it in the following manner ^{[11]}:

**(**)

where

**(**)

**(**)

**(**)

**(**)

and

**(**)

In matrix form, equation (**C-68**) now is written as follows:

**(**)

Consider an initial estimate of the parameter vector **p** : (…, *s _{wj}*, …,

*s*, …,

_{wi}*s*). We want to minimize the difference between the observed and the modeled times by iteratively perturbing the initial estimate of the parameter vector. A change Δ

_{b}*p*in the paramater vector will change the modeled times as follows:

**(**)

The error in modeling the traveltimes is given by

**(**)

Substitute equation (**C-71**) into equation (**C-72**) to obtain

**(**)

Now define the difference Δ*t _{ij}* between the observed traveltimes

*t*and the initial estimate of the modeled traveltimes and rewrite equation (

_{ij}**C-73**) to get

**(**)

The second term on the right is the amount of change in Δ*t _{ij}* as a result of the change in the parameter Δ

*p*. Define this term as and rewrite equation (

**C-74**) once more to obtain

**(**)

where

**(**)

The derivatives in equation (**C-76**) can be computed by differentiating the model equation (**C-68**) with respect to each of the parameters:

**(**)

**(**)

and

**(**)

These then are substituted back into the right-hand side of equation (**C-76**) to get

**(**)

Examine the structure of equation (**C-78**) and note that, instead of modeling refraction traveltimes by way of equation (**C-68**), we can model the change in traveltimes by way of equation (**C-78**), and thus estimate the near-surface parameters. Hence, for *n* shot/receiver stations, the parameter vector is **Δp** : (Δ*s _{w}*

_{1}, Δ

*s*

_{w}_{2}, …, Δ

*s*, Δ

_{wn}*s*). We refer to the scheme based on equation (

_{b}**C-78**) as the

*iterative*GLI solution. For

*m*picks of

*t*, and

_{ij}*n*+ 1 parameters

**Δp**: (Δ

*s*

_{w}_{1}, Δ

*s*

_{w}_{2}, …, Δ

*s*, Δ

_{wn}*s*), we have the following set of equations:

_{b}

**(**)

We write these equations in matrix notation as in equation (**C-59**) to obtain

**(**)

or

**(**)

where **Δt′** is the column vector of *m*-length on the left-hand side of equation (**C-79**), **L** is the sparse matrix in the same equation with nonzero elements as the partial derivatives and with dimensions *m* × (*n* + 1), and **p** is the column vector of (*n* + 1)-length on the right-hand side of equation (**C-79**). Except for the three elements in each row, the **L** matrix contains zeros.

The GLI solution to equation (**C-81**) satisfies the requirement that the energy of the error vector

**(**)

is minimum and is given by

**(**)

where **Δt** denotes the column vector of *m*-length that represents the difference between the observed (picked) refracted arrival times and the initial estimate of the modeled times, and *T* denotes matrix transposition.

In summary, the iterative scheme for refraction statics is as follows:

- Specify a flat datum to which shot and receivers are to be lowered along vertical raypaths.
- Also specify a set of initial model parameters
**p**: (*s*_{w}_{1},*s*_{w}_{2}, …,*s*,_{wn}*s*)._{b} - Compute Δ
*t*, the time difference between the picked (observed) times_{ij}*t*and the initial modeled times ._{ij} - Estimate the change in parameters:
**Δp**: (Δ*s*_{w}_{1}, Δ*s*_{w}_{2}, …, Δ*s*, Δ_{wn}*s*), by way of the GLI solution given by equation (_{b}**C-83**). - Update the paramater vector
**p**+**Δp**, and compute new modeled times . - Iterate steps (c), (d), and (e) to get a final estimate of the parameter vector
**p**.

The parameter vector for the iterative scheme described above comprises laterally varying weathering velocity and a constant bedrock velocity. Depth to the refractor is assumed to be known. An alternative parameterization of the near-surface model may involve estimating a laterally varying weathering velocity and depth to the refractor while assuming a known value for the bedrock velocity.

Rearrange the terms in equation (**C-62**):

**(**)

so that

**(**)

Equation (**C-84**) implies that, in this scheme, we deal with picked times that have been linear-moveout corrected using the assumed value for the bedrock velocity. By way of equations (**C-55**) and (**C-56**), equation (**C-85**) is expressed as follows:

**(**)

Compute the partial derivatives of the parameters to be estimated — depth to refractor and weathering velocity at shot and receiver locations:

**(**)

**(**)

**(**)

and

**(**)

Now rewrite equation (**C-76**) in terms of the new variable *t"* to obtain

**(**)

and substitute equations (2-49a, b, c, d) to get the model equation

**(**)

Finally, similar to equation (**C-79**), we have the following set of equations:

**(**)

Follow the steps which involve equations (**C-80**) through (**C-83**) to estimate the parameter vector **p** : (…, *z _{j}*, …,

*z*, …,

_{i}*s*, …,

_{wj}*s*).

_{wi}## C.10 *L*_{1}-norm refraction statics

The generalized linear inversion method applied to residual and refraction statics is based on minimization of the quantity

**(**)

where *t _{ij}* are the actual traveltime picks and are the modeled traveltimes defined by equation (

**C-24a**) for residual statics and (

**C-54**) for refraction statics.

The minimization norm defined by equation (**C-91**) is formally referred to as the *L*_{2} norm. For statics applications, it may be desirable to use the *L*_{1} minimization norm defined by

**(**)

Outliers in picked times can cause biased results from the *L*_{2}-norm, whereas, they may be better handled by the *L*_{1}-norm schemes. The *L*_{2}-norm solutions to equations (**C-24a**) and (**C-54**) are expressed by the generalized linear inversion formula (equations **C-28** or **C-61**). Whereas, the *L*_{1}-norm solution is based on linear programming techniques ^{[15]} and will not be dealt with here. Instead, we shall refer to model experiments and a field data example for the *L*_{1}-norm refraction statics corrections.

**Figure C-6**CMP stack associated with the model data in Figure 3.4-14 — (a) with refraction statics corrections using the*L*_{1}-norm, (b) section as in (a) after residual statics corrections using the*L*_{1}-norm. Compare with the*L*_{2}-norm results in Figure 3.4-17.**Figure C-7**Summary of the*L*_{1}-norm solution for refraction statics associated with the CMP stacked section in Figure 3.4-15a. Plot direction is the same as that in Figure 3.4-15. Except in frame 1, shot attributes are denoted with × and receiver attributes are denoted with vertical bars. Estimated parameters are intercept time anomalies and are plotted in frame 1 with no distinction made between shot and receiver locations.**Figure C-8**CMP stack associated with the model data in Figure 3.4-21 — (a) with refraction statics corrections using the*L*_{1}-norm, (b) section as in (a) after residual statics corrections using the*L*_{1}-norm. Compare with the*L*_{2}-norm results in Figure 3.4-26.**Figure C-9**Summary of the*L*_{1}-norm solution for refraction statics associated with the CMP stacked section in Figure 3.4-26a. Plot direction is the same as that in Figure 3.4-26. Except in frame 1, shot attributes are denoted with × and receiver attributes are denoted with vertical bars. Estimated parameters are intercept time anomalies and are plotted in frame 1 with no distinction made between shot and receiver locations.**Figure C-10**The CMP stacked section of the data in Figure 3.4-29 (a) after refraction statics using an*L*_{1}-norm scheme for refraction statics, and (b) followed by residual statics corrections. Compare with Figure 3.4-29 and note the significant improvement of CMP stacking as a result of resolving both long- and short-wavelength statics anomalies by refraction and residual statics corrections, respectively. Also, compare with Figure 3.4-32 and note that the results of*L*_{1}- and*L*_{2}-norm statics solutions, in this case, yield very similar results.

Consider the single-layered near-surface model shown in Figure 3.4-14. Figure C-6 shows results of refraction and residual statics corrections using the *L*_{1}-norm. Note that these results are comparable to those obtained from *L*_{2}-norm minimization (Figure 3.4-17).

The results of the *L*_{1}-norm statics estimates are summarized in Figure C-7. Compare with those of the *L*_{2}-norm solution in Figure 3.4-18 and note that the statics solutions exhibit minor differences, although the resulting stacked sections are almost identical. Frame 1 shows the estimated intercept time anomalies (equations **53a**, **53b**), as a function of the shot-receiver station number. Frame 2 shows the pick fold, namely the number of picks in each shot (denoted by ×) and receiver (denoted by the vertical bars) gather.

A quantitative measure of the accuracy of the *L*_{1}-norm solution to refraction statics is the sum of the differences between the observed picks *t _{ij}* and the modeled traveltimes (equation

**C-91b**) over each shot and receiver gather. These residual time differences are plotted in frame 3 of Figure C-7. Large residuals often are related to bad picks. Nevertheless, even with good picks, there may be large residuals attributable to inappropriateness of the model assumed for the near-surface.

Figure C-7 also shows the estimated weathering thicknesses at all shot-receiver stations (frame 4). Finally, the computed statics and the near-surface model are shown in frames 5 and 6, respectively.

Now consider the multilayered near-surface model shown in Figure 3.4-21. Figure C-8 shows results of refraction and residual statics corrections using the *L*_{1}-norm solution. Again, results are comparable to those obtained from *L*_{2}-norm minimization (Figure 3.4-26).

Results of the *L*_{1}-norm refraction statics solution are summarized in Figure C-9. Compare with those of the *L*_{2}-norm solution in Figure 3.4-27 and note that the statics solutions and the resulting stacked sections are virtually identical. Although the actual near-surface model consists of several layers (Figure 3.4-21a), the *L*_{1}-norm solution is based on a single layer with a constant velocity of 1400 m/s as for the *L*_{2}-norm solution (Figure 3.4-27). In Figure C-9, Frame 1 shows the estimated intercept time anomalies (equation **53a**, **C-91b**) over each shot and receiver gather is shown in frame 3. Large residuals, in this case, are attributable to inappropriateness of the model assumed for the near-surface. Figure C-9 also shows the estimated weathering thicknesses at all shot-receiver stations (frame 4). Finally, the computed statics and the nearsurface model are shown in frames 5 and 6, respectively. Whatever the inversion norm, the important point is how the near-surface model is parameterized and how close it is to the real situation.

Figure C-10 shows the results of refraction and residual statics corrections applied to field data as in Figure 3.4-29 using the *L*_{1}-norm solution. The results are comparable to those obtained from *L*_{2}-norm minimization (Figure 3.4-32).

## Figures

**Figure 3.4-11**(a) Geometry for refracted arrivals. Here,*v*= weathering velocity,_{w}*v*= bedrock velocity,_{b}*z*= depth to the refractor equivalent to the base of the weathering layer,_{w}*θ*= critical angle, and_{c}*x*= crossover distance. The direct wave arrival has a slope equal to 1/_{c}*v*and the refracted wave arrival has a slope equal to 1/_{w}*v*. (b) A shot record that exhibits the direct wave and the refracted wave depicted in (a). (c) Geometry for a dipping refractor with forward traveltime profile associated with the direct wave and refracted wave arrivals, and (d) with both forward and reverse traveltime profiles. See text for details._{b}**Figure 3.4-12**(a) Geometry for the plus-minus method. (b) Geometry for the generalized reciprocal method. Here,*z*is the depth to the refractor at the surface station where the plus-minus times as for (a) and intercept times as for (b) are to be estimated,_{w}*v*is the weathering velocity, and_{w}*θ*is the critical angle of refraction._{c}**Figure 3.4-13**Geometry of refracted arrival used in deriving the least-squares solution for intercept times. Here,*S*and_{j}*R*are source and receiver stations, respectively;_{i}*θ*is the critical angle of refraction,_{c}*z*and_{j}*z*are depths to the bedrock at source and receiver locations, and_{i}*v*and_{w}*v*are weathering and bedrock velocities, respectively._{b}**Figure 3.4-17**CMP stack associated with the model data in Figure 3.4-14 — (a) with refraction statics corrections using the GLI solution, (b) section as in (a) after residual statics corrections. Datum level is 0 m in both sections, and the weathering velocity is assumed to be 1200 m/s.**Figure 3.4-18**Summary of the variable-thickness GLI solution for refraction statics associated with the CMP stacked section in Figure 3.4-17a. See text for details. Plot direction is the same as that in Figure 3.4-17. Except in frame 1, shot attributes are denoted with × and receiver attributes are denoted with vertical bars. Estimated parameters from equation (**C-54**) are plotted in frame 1 with no distinction made between shot and receiver locations.**Figure 3.4-26**CMP stack associated with the model data in Figure 3.4-21 — (a) with refraction statics corrections using the GLI solution, (b) section as in (a) after residual statics corrections. Datum level is 0 m in both sections, and the weathering velocity is assumed to be 1400 m/s.**Figure 3.4-27**Summary of the variable-thickness GLI solution for refraction statics associated with the CMP stacked section in Figure 3.4-26a. Plot direction is the same as that in Figure 3.4-26. Except in frame 1, shot attributes are denoted with × and receiver attributes are denoted with vertical bars. Estimated parameters from equation (**C-54**) are plotted in frame 1 with no distinction made between shot and receiver locations.**Figure 3.4-29**(a) A CMP stack with elevation statics corrections, (b) same as in (a) with residual statics corrections. (Data courtesy Nederlandse Aardolie Maatschappij B.V.)**Figure 3.4-32**(a) A CMP stack with refraction statics applied using the least-squares method (compare with Figures 3.4-29a and Figure 3.4-31a), (b) same as in (a) with residual statics corrections.

## Equations

**(**)

**(**)

**(**)

**(**)

**(**)

**(**)

**(**)

**(**)

**(**)

**(**)

**(**)

**(**)

## See also

- Introduction to velocity analysis and statics corrections
- Normal moveout
- Velocity analysis
- Residual statics corrections
- Refraction statics corrections
- Exercises

## References

- ↑ Taner and Koehler, 1969, Taner, M. T. and Koehler, F., 1969, Velocity spectra — digital computer derivation and applications of velocity functions: Geophysics, 39, 859–881.
- ↑
^{2.0}^{2.1}^{2.2}^{2.3}Castle (1994), Castle, R. J., 1994, A theory of normal moveout: Geophysics, 59, 983–999. - ↑
^{3.0}^{3.1}Slotnick, 1959, Slotnick, M. M., 1959, Lessons in seismic computing: Soc. Expl. Geophys. - ↑
^{4.0}^{4.1}Grant and West, 1965 - ↑
^{5.0}^{5.1}De Bazelaire (1988), De Bazelaire, E., 1988, Normal moveout revisited: Inhomogeneous media and curved interfaces: Geophysics, 53, 143–157. - ↑ Taner et al., 1974, Taner, M. T., Koehler, F., and Alhilali, K. A., 1974, Estimation and correction of near-surface time anomalies: Geophysics, 41, 441–463.
- ↑ Wiggins et al., 1976, Wiggins, R. A., Larner, K. L., and Wisecup, R. D., 1976, Residual statics analysis as a general linear inverse problem: Geophysics, 41, 922–938.
- ↑
^{8.0}^{8.1}Lines and Treitel, 1984, Lines, L. R. and Treitel, S., 1984, Tutorial: A review of least-squares inversion and its application to geophyiscal problems: Geophys. Prosp., 32, 159–186. Cite error: Invalid`<ref>`

tag; name "ch03r21" defined multiple times with different content - ↑
^{9.0}^{9.1}Hampson and Russell, 1984, Hampson, D. and Russell, B., 1984, First-break interpretation using generalized inversion: J. Can. Soc. Explor. Geophys., 20, 40–54. - ↑
^{10.0}^{10.1}Schneider and Kuo, 1985, Schneider, W. A. and Kuo, S., 1985, Refraction modeling for statics corrections: 55th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 295–299. - ↑
^{11.0}^{11.1}^{11.2}De Amorim et al., 1987, De Amorim, W. N., Hubral, P. and Tygel, M., 1987, Computing field statics with the help of seismic tomography: Geophys. Prosp., 35, 907–919. - ↑ Baixas and Du Pont, 1988, Baixas, F. and Du Pont, R., 1988, Practical view of 3-D refraction statics: 58th Ann. Internat. Mtg., Soc. Expl. Geophs., Expanded Abstracts, 787–790.
- ↑ Kircheimer, 1988, Kircheimer, F., 1988, 3-D Refraction statics by weighted least-squares inversion: 59th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 794–797.
- ↑ Farrell and Euwema, 1984, Farrell, R. C. and Euwema, R. N., 1984, Refraction Statics: Proc. Inst. Electr. Electron. Eng., 72, 1316–1329.
- ↑ Press et al., 1987, Press, W. H., B. P. Flannery, S. A. Teukolsky and W. A. Vettering, 1987: Numerical recipes — The art of scientific computing, Cambridge University Press.