# Topics in moveout and statics corrections

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

## 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]:

 ${\displaystyle t^{2}=C_{0}+C_{1}x^{2}+C_{2}x^{4}+C_{3}x^{6}+\cdots ,}$ (C-1)

where x is the offset, and

 ${\displaystyle C_{0}=t_{0}^{2},}$ (C-2a)

where t0 is the two-way zero-offset traveltime

 ${\displaystyle C_{1}={\frac {1}{\mu _{2}}},}$ (C-2b)

 ${\displaystyle C_{2}={\frac {1}{4}}{\frac {\mu _{2}^{2}-\mu _{4}}{t_{0}^{2}\mu _{2}^{4}}},}$ (C-2c)

and

 ${\displaystyle C_{3}={\frac {2\mu _{4}^{2}-\mu _{2}\mu _{6}-\mu _{2}^{2}\mu _{4}}{t_{0}^{4}\mu _{2}^{7}}},}$ (C-2d)

with the additional definitions of the terms [2]

 ${\displaystyle \mu _{j}={\frac {1}{t_{0}}}\sum \limits _{i=1}^{N}v_{i}^{j}\Delta \tau _{i}.}$ (C-3a)

Note that

 ${\displaystyle \mu _{2}=v_{rms}^{2}}$ (C-3b)

since

 ${\displaystyle v_{rms}^{2}={\frac {1}{t_{0}}}\sum _{i=1}^{N}v_{i}^{2}\Delta \tau _{i},}$ (C-3c)

where Δti is the vertical two-way traveltime through the i-th layer, v is the velocity, and ${\displaystyle {t_{0}}=\sum \nolimits _{i=1}^{N}\Delta \tau _{i}.}$ 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

 ${\displaystyle t^{2}=C_{0}+C_{1}x^{2}+C_{2}x^{4}.}$ (C-4)

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

 ${\displaystyle t^{2}=t_{0}^{2}+{\frac {x^{2}}{v_{rms}^{2}}}+C_{2}x^{4}.}$ (C-5)

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

 ${\displaystyle t^{2}=t_{0}^{2}+{\frac {x^{2}}{v_{rms}^{2}}}.}$ (C-6)

Compute the velocity spectrum using equation (C-6), and pick an initial velocity function vrms(t0). Then, use this picked velocity function in equation (C-5) and compute a velocity spectrum for the parameter C2. Pick a function C2(t0) and use in equation (C-5) to recompute the velocity spectrum. Finally, pick an updated velocity function vrms(t0) 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]

 ${\displaystyle t=t_{0}\left(1-{\frac {1}{S}}\right)+{\sqrt {\left({\frac {t_{0}}{S}}\right)^{2}+{\frac {x^{2}}{Sv_{rms}^{2}}}}},}$ (C-7)

where S is a constant of the form

 ${\displaystyle S={\frac {\mu _{4}}{\mu _{2}^{2}}}.}$ (C-8)

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 = t0(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:

1. The traveltime trajectory is symmetric with respect to the time axis.
2. The instantaneous velocity dx/dt is never zero.
3. 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 vrms(t0). Then, use this picked velocity function in equation (C-7) and compute a velocity spectrum for the parameter S. Pick a function S(t0) and use in equation (C-7) to recompute the velocity spectrum. Finally, pick an updated velocity function vrms(t0) from this velocity spectrum.

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

 ${\displaystyle t=(t_{0}-t_{p})+{\sqrt {t_{p}^{2}+{\frac {x^{2}}{v_{s}^{2}}}}},}$ (C-9)

where t0 is the two-way zero-offset time, tp is related to the time at which the asymptotes of the hyperbolic traveltime trajectory converge (Figure C-2), and vs is the reference velocity assigned to the layer below the recording surface (not the near-surface layer). When tp = t0, equation (C-9) reduces to the small-spread hyperbolic equation (C-6).

Figure C-1  Traveltime trajectories based on (top) the hyperbolic equation (C-23) and (bottom) the time-shifted hyperbolic equation (C-14). Compare with the true traveltime trajectory associated with a layered model [2].

## 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

 ${\displaystyle t^{2}=t_{0}^{2}+{\frac {x^{2}}{v^{2}}},}$ (C-10)

where t is the two-way traveltime associated with a source-receiver separation x, t0 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 T0 = 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

 ${\displaystyle (t+T)^{2}=(t_{0}+T+\Delta T)^{2}+{\frac {x^{2}}{v^{2}}}.}$ (C-11)
Figure C-2  The traveltime trajectory associated with the moveout equation (5c) [5].

Expand the terms on both sides:

 ${\displaystyle t^{2}+2tT+T^{2}=t_{0}^{2}+2t_{0}(T+\Delta T)+(T+\Delta T)^{2}+{\frac {x^{2}}{v^{2}}}.}$ (C-12a)

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

 ${\displaystyle 2tT+T^{2}=2t_{0}(T+\Delta T)+(T+\Delta T)^{2}.}$ (C-12b)

Simplify and rearrange the terms

 ${\displaystyle 2(t-t_{0})T=2(t_{0}+T)\Delta T+\Delta T^{2}.}$ (C-12c)

Now, ignore the second term on the right-hand side of the equation and observe that ΔtNMO = tt0 to obtain

 ${\displaystyle \Delta t_{NMO}T=(t_{0}+T)\Delta T.}$ (C-12d)

Assume that t0 >> T and rearrange the terms to obtain a relationship for change in the period of the wavelet as a result of moveout correction:

 ${\displaystyle {\frac {\Delta T}{T}}={\frac {\Delta t_{NMO}}{t_{0}}}.}$ (C-13)

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

 ${\displaystyle T={\frac {1}{f}},}$ (C-14)

and obtain

 ${\displaystyle \Delta T=-{\frac {1}{f^{2}}}\Delta f.}$ (C-15)

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

 ${\displaystyle {\frac {\Delta f}{f}}={\frac {\Delta t_{NMO}}{t_{0}}}.}$ (C-16)

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 : (xF, zF) and G : (xG, zG), so that

 ${\displaystyle (FG)^{2}=(x_{F}-x_{G})^{2}+(z_{F}-z_{G})^{2}.}$ (C-17)

From the geometry of Figure C-3, the coordinates of F : (xF, zF) are

${\displaystyle x_{F}=OF\cos 2\phi }$

and

${\displaystyle z_{F}=OF\sin 2\phi ,}$

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

 ${\displaystyle x_{F}=\left(OM+{\frac {x}{2}}\right)\cos 2\phi }$ (C-18a)

and

 ${\displaystyle z_{F}=\left(OM+{\frac {x}{2}}\right)\sin 2\phi .}$ (C-18b)

Again, from the geometry of Figure C-3, the coordinates of G : (xG, zG) are

 ${\displaystyle x_{G}=OM-{\frac {x}{2}},}$ (C-19a)

and

 ${\displaystyle z_{G}=0.}$ (C-19b)

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

 ${\displaystyle (FG)^{2}=\left[(OM+{\frac {x}{2}})\cos 2\phi -(OM-{\frac {x}{2}})\right]^{2}+\left[(OM+{\frac {x}{2}})\sin 2\phi \right]^{2},}$ (C-20a)

then use the trigonometric relation cos 2ϕ = 2 cos2 ϕ − 1 and apply some algebra to obtain

 ${\displaystyle (FG)^{2}=\left(OM+{\frac {x}{2}}\right)^{2}+\left(OM-{\frac {x}{2}}\right)^{2}-2\left[(OM)^{2}-{\frac {x^{2}}{4}}\right](2\cos ^{2}\phi -1).}$ (C-20b)

Further algebraic manipulation yields

 ${\displaystyle (FG)^{2}=4(OM)^{2}\sin ^{2}\phi +x^{2}\cos ^{2}\phi .}$ (C-21)

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

 ${\displaystyle MN=OM\sin \phi .}$ (C-22)

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

 ${\displaystyle t^{2}=t_{0}^{2}+{\frac {x^{2}\cos ^{2}\phi }{v^{2}}},}$ (C-23a)

and the moveout velocity

 ${\displaystyle v_{NMO}={\frac {v}{\cos \phi }}.}$ (C-23b)

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

Figure C-3  Geometry of a dipping reflector used in deriving the equations in Section C.3.

## C.4 Traveltime decomposition for residual statics estimation

We want to model traveltime deviations ${\displaystyle t'_{ij}}$ associated with a reflection event on moveout-corrected CMP gathers by the following equation [6]; [7]:

 ${\displaystyle t'_{ij}=s_{j}+r_{i}+G_{k}+M_{k}x_{ij}^{2},}$ (C-24a)

where sj is the residual statics shift at the jth source location, ri is the residual statics shift at the ith receiver location, Gk is the structure term at the kth midpoint location, [k = (i + j)/2], ${\displaystyle {M_{k}}{x_{ij}^{2}}}$ is the residual moveout at the kth midpoint location.

For m picks of tij, and ns shot locations, nr receiver locations, nG midpoint locations, we have the following set of equations:

 ${\displaystyle {\begin{pmatrix}\vdots \\\vdots \\\vdots \\t'_{ij}\\\vdots \\\vdots \\\vdots \\\end{pmatrix}}={\begin{pmatrix}\ \\\ \\\cdots 1\ \cdots 1\ \cdots 1\ \cdots x_{ij}^{2}\cdots \\\ \\\ \\\end{pmatrix}}{\begin{pmatrix}\vdots \\s_{j}\\\vdots \\r_{i}\\\vdots \\G_{k}\\\vdots \\M_{k}\\\vdots \\\end{pmatrix}}}$ (C-24b)

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

 ${\displaystyle t'_{ij}=s_{j}+r_{i}+Mx_{ij}^{2}.}$ (C-25a)

Equation (C-24b) is modified, accordingly

 ${\displaystyle {\begin{pmatrix}\vdots \\\vdots \\\vdots \\t'_{ij}\\\vdots \\\vdots \\\vdots \\\end{pmatrix}}={\begin{pmatrix}\ \\\cdots \quad 1\quad \cdots \quad 1\quad \cdots \quad x_{ij}^{2}\\\ \\\end{pmatrix}}{\begin{pmatrix}\vdots \\s_{j}\\\vdots \\s_{i}\\\vdots \\M\\\end{pmatrix}}}$ (C-25b)

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

 ${\displaystyle \mathbf {t'} =\mathbf {Lp} ,}$ (C-26)

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 × (ns + nr + nG + nG) in case of equation (C-24b) and m × (ns + nr + 1) in case of equation (C-25b), and p is the column vector of (ns + nr + nG + nG)-length in case of equation (C-24b) and (ns + nr + 1)-length in case of equation (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

 ${\displaystyle \mathbf {e} =\mathbf {t-t'} }$ (C-27a)

is minimum [8]. The cumulative error energy is

 ${\displaystyle E=\mathbf {e^{T}e} .}$ (C-27b)

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

 ${\displaystyle E=(\mathbf {t} -\mathbf {Lp} )^{T}(\mathbf {t} -\mathbf {Lp} ).}$ (C-27c)

Minimization of E with respect to p requires that

 ${\displaystyle {\frac {\partial E}{\partial s_{j}}}={\frac {\partial E}{\partial r_{i}}}={\frac {\partial E}{\partial G_{k}}}={\frac {\partial E}{\partial M_{k}}}=0,}$ (C-27d)

which, in the case of equation (C-24b), yields ns + nr + nG + nG equations and that many unknowns. These equations can be solved for the residual statics associated with ns source locations, nr receiver locations, nG structural terms, and nG residual moveout terms.

Finally, the solution that satisfies the minimization requirement given by equation (C-27d) follows [8]

 ${\displaystyle \mathbf {p} =(\mathbf {L^{T}L} )^{-1}\mathbf {L^{T}t} ,}$ (C-28)

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

 ${\displaystyle s_{j}^{m}={\frac {1}{n_{r}}}\sum _{i}^{n_{r}}\{t_{ij}-G_{k}^{m-1}-M_{k}^{m-1}x_{ij}^{2}-r_{i}^{m-1}\},}$ (C-29a)

 ${\displaystyle r_{i}^{m}={\frac {1}{n_{s}}}\sum _{j}^{n_{s}},\{t_{ij}-G_{k}^{m-1}-M_{k}^{m-1}x_{ij}^{2}-s_{j}^{m-1}\},}$ (C-29b)

 ${\displaystyle G_{k}^{m}={\frac {1}{n_{h}}}\sum _{l}^{n_{h}},\{t_{ij}-s_{j}^{m-1}-M_{k}^{m-1}x_{ij}^{2}-r_{i}^{m-1}\},}$ (C-29c)

and

 ${\displaystyle M_{k}^{m}={\frac {1}{n_{G}}}\sum _{k}^{n_{G}}{\frac {1}{x_{ij}^{2}}}\{t_{ij}-s_{j}^{m-1}-G_{k}^{m-1}-r_{i}^{m-1}\},}$ (C-29d)

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 nh 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

 ${\displaystyle s_{j}^{m}={\frac {1}{n_{r}}}\sum _{i}^{n_{r}}\{t_{ij}\}-{\frac {1}{n_{r}}}\sum _{i}^{n_{r}}\{G_{k}^{m-1}-M_{k}^{m-1}x_{ij}^{2}-r_{i}^{m-1}\},}$ (C-30a)

 ${\displaystyle r_{i}^{m}={\frac {1}{n_{s}}}\sum _{j}^{n_{s}}\{t_{ij}\}-{\frac {1}{n_{s}}}\sum _{j}^{n_{s}}\{G_{k}^{m-1}-M_{k}^{m-1}x_{ij}^{2}-s_{j}^{m-1}\},}$ (C-30b)

 ${\displaystyle G_{k}^{m}={\frac {1}{n_{h}}}\sum _{l}^{n_{h}}\{t_{ij}\}-{\frac {1}{n_{h}}}\sum _{l}^{n_{h}}\{s_{j}^{m-1}-M_{k}^{m-1}x_{ij}^{2}-r_{i}^{m-1}\},}$ (C-30c)

and

 ${\displaystyle M_{k}^{m}={\frac {1}{n_{G}}}\sum _{k}^{n_{G}}{\frac {1}{x_{ij}^{2}}}\{t_{ij}\}-{\frac {1}{n_{G}}}\sum _{k}^{n_{G}}{\frac {1}{x_{ij}^{2}}}\{s_{j}^{m-1}-G_{k}^{m-1}-r_{i}^{m-1}\}.}$ (C-30d)

This modification enables us to compute and store the sum of the traveltime deviations ∑tij 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 sj, ri, Gk, and Mk from the mth 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 zw (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

 ${\displaystyle t={\frac {SB}{v_{w}}}+{\frac {BC}{v_{b}}}+{\frac {CR}{v_{w}}}.}$ (C-31a)

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

 ${\displaystyle t={\frac {z_{w}}{v_{w}\cos \theta _{c}}}+{\frac {x-2z_{w}\tan \theta _{c}}{v_{b}}}+{\frac {z_{w}}{v_{w}\cos \theta _{c}}},}$ (C-31b)

where x is the shot-receiver separation, zw is the depth to bedrock, and vw and vb are weathering and bedrock velocities, respectively. The critical angle of refraction θc is related to the velocities by the relation

 ${\displaystyle \sin \theta _{c}={\frac {v_{w}}{v_{b}}}.}$ (C-32)

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

 ${\displaystyle t={\frac {2z_{w}}{v_{w}\cos \theta _{c}}}-{\frac {2z_{w}\sin \theta _{c}}{v_{b}\cos \theta _{c}}}+{\frac {x}{v_{b}}}.}$ (C-33a)

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

 ${\displaystyle t={\frac {2z_{w}{\sqrt {v_{b}^{2}-v_{w}^{2}}}}{v_{b}v_{w}}}+{\frac {x}{v_{b}}}.}$ (C-33b)

Note that this is the equation of a line

 ${\displaystyle t=t_{i}+{\frac {x}{v_{b}}},}$ (C-34)

with its slope given by 1/vb — inverse of the bedrock velocity, and the intercept time ti given by

 ${\displaystyle t_{i}={\frac {2z_{w}{\sqrt {v_{b}^{2}-v_{w}^{2}}}}{v_{b}v_{w}}}.}$ (C-35)

By measuring the slope of the line associated with refracted arrivals in Figure 3.4-11a, we estimate the bedrock velocity vb. We also estimate the intercept time ti by extending the line to x = 0. Finally, we assume a value for the weathering velocity vw. Then, equation (C-35) can be rearranged to compute the depth to bedrock zw as

 ${\displaystyle z_{w}={\frac {v_{b}v_{w}t_{i}}{2{\sqrt {v_{b}^{2}-v_{w}^{2}}}}}.}$ (C-36)

This is equation (41a) in the text.

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

 ${\displaystyle t={\frac {x}{v_{w}}},}$ (C-37)

whereas equation (C-34) describes the refracted arrivals. At the crossover distance xc, we equate the arrival times in equations (C-34) and (C-37) as

 ${\displaystyle t_{i}+{\frac {x_{c}}{v_{b}}}={\frac {x_{c}}{v_{w}}}.}$ (C-38a)

By substituting equation (C-35), we obtain

 ${\displaystyle z_{w}={\frac {v_{b}v_{w}x_{c}}{2{\sqrt {v_{b}^{2}-v_{w}^{2}}}}}\left({\frac {1}{v_{w}}}-{\frac {1}{v_{b}}}\right).}$ (C-38b)

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

 ${\displaystyle z_{w}={\frac {1}{2}}{\sqrt {\frac {v_{b}-v_{w}}{v_{b}+v_{w}}}}x_{c}.}$ (C-39)

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

 ${\displaystyle t^{-}={\frac {SA}{v_{w}}}+{\frac {AB}{v_{b}}}+{\frac {BR}{v_{w}}}.}$ (C-40a)

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

 ${\displaystyle t^{-}={\frac {z_{wS}}{v_{w}\cos \theta _{c}}}+{\frac {x\cos \varphi -z_{wS}\tan \theta _{c}-(z_{wS}+x\sin \varphi )\tan \theta _{c}}{v_{b}}}+{\frac {z_{wS}+x\sin \varphi }{v_{w}\cos \theta _{c}}}.}$ (C-40b)

Apply some algebra to obtain

 ${\displaystyle t^{-}={\frac {2z_{wS}\cos \theta _{c}\cos \varphi }{v_{w}}}+{\frac {x\sin(\theta _{c}+\varphi )}{v_{w}}}.}$ (C-41a)

This is the equation of a line

 ${\displaystyle t^{-}=t_{i}^{-}+{\frac {x}{v_{b}^{-}}},}$ (C-41b)

with its inverse slope defined as

 ${\displaystyle v_{b}^{-}={\frac {v_{w}}{\sin(\theta _{c}+\varphi )}}}$ (C-41c)

and the intercept time given by

 ${\displaystyle t_{i}^{-}={\frac {2z_{wS}\cos \theta _{c}\cos \varphi }{v_{w}}}.}$ (C-41d)

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

 ${\displaystyle t^{+}={\frac {2z_{wR}\cos \theta _{c}\cos \varphi }{v_{w}}}+{\frac {x\sin(\theta _{c}-\varphi )}{v_{w}}}.}$ (C-42a)

Again, this is the equation of a line

 ${\displaystyle t^{+}=t_{i}^{+}+{\frac {x}{v_{b}^{+}}},}$ (C-42b)

with its inverse slope defined as

 ${\displaystyle v_{b}^{+}={\frac {v_{w}}{\sin(\theta _{c}-\varphi )}},}$ (C-42c)

and the intercept time given by

 ${\displaystyle t_{i}^{+}={\frac {2z_{wR}\cos \theta _{c}\cos \varphi }{v_{w}}}.}$ (C-42d)

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

 ${\displaystyle \theta _{c}+\varphi =\sin ^{-1}{\frac {v_{w}}{v_{b}^{-}}}}$ (C-43a)

and

 ${\displaystyle \theta _{c}-\varphi =\sin ^{-1}{\frac {v_{w}}{v_{b}^{+}}}.}$ (C-43b)

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

 ${\displaystyle \varphi ={\frac {1}{2}}\left[\sin ^{-1}{\frac {v_{w}}{v_{b}^{-}}}-\sin ^{-1}{\frac {v_{w}}{v_{b}^{+}}}\right].}$ (C-44)

This is equation (46a) of the text. Measure the slopes of the downdip and updip refracted arrivals in Figure 3.4-11d — ${\displaystyle {v_{b}^{-}}\ {\text{and }}{v_{b}^{+}},}$ respectively. Then, assume a value for the weathering velocity vw. By direct substitution into equation (C-44), compute the refractor dip φ.

To obtain the bedrock velocity vb, first, rewrite slopes from equations (C-41c) and (C-42c) as

 ${\displaystyle {\frac {1}{v_{b}^{-}}}={\frac {\sin \theta _{c}\cos \varphi +\sin \varphi \cos \theta _{c}}{v_{w}}}}$ (C-45a)

and

 ${\displaystyle {\frac {1}{v_{b}^{+}}}={\frac {\sin \theta _{c}\cos \varphi -\sin \varphi \cos \theta _{c}}{v_{w}}}.}$ (C-45b)

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

 ${\displaystyle v_{b}={\frac {2\cos \varphi }{\left({\frac {1}{v_{b}^{-}}}+{\frac {1}{v_{b}^{+}}}\right)}}.}$ (C-46)

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 vb from equation (C-46), and the measurements for the intercept times ${\displaystyle {t_{i}^{-}}\ {\text{and }}{t_{i}^{+}}}$ from the arrival times in Figure 3.4-11d into equations (C-41d) and (C-42d) to obtain

 ${\displaystyle z_{wS}={\frac {v_{b}v_{w}t_{i}^{-}}{2\cos \varphi {\sqrt {v_{b}^{2}-v_{w}^{2}}}}}}$ (C-47a)

and

 ${\displaystyle z_{wR}={\frac {v_{b}v_{w}t_{i}^{+}}{2\cos \varphi {\sqrt {v_{b}^{2}-v_{w}^{2}}}}}.}$ (C-47b)

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

 ${\displaystyle t_{+}=t_{ABCD}+t_{DEFG}-t_{ABFG}}$ (C-48a)

and

 ${\displaystyle t_{-}=t_{ABCD}-t_{DEFG}+t_{ABFG}.}$ (C-48b)

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

 ${\displaystyle t_{+}=2\left({\frac {CD}{v_{w}}}-{\frac {CH}{v_{b}}}\right).}$ (C-49a)

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 zw as

 ${\displaystyle t_{+}=2\left({\frac {z_{w}}{v_{w}\cos \theta }}-{\frac {z_{w}\tan \theta }{v_{b}}}\right).}$ (C-49b)

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 zw, weathering velocity vw, and bedrock velocity vb:

 ${\displaystyle t_{+}={\frac {2z_{w}{\sqrt {v_{b}^{2}-v_{w}^{2}}}}{v_{w}v_{b}}}.}$ (C-48c)

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

 ${\displaystyle t_{-}={\frac {2CD}{v_{w}}}+{\frac {2BC}{v_{b}}}+{\frac {CE}{v_{b}}}.}$ (C-49a)

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

 ${\displaystyle t_{-}={\frac {2z_{w}}{v_{w}\cos \theta }}+{\frac {2BC}{v_{b}}}+{\frac {2z_{w}\tan \theta }{v_{b}}}.}$ (C-49b)

By some algebraic manipulation, we get

 ${\displaystyle t_{-}={\frac {2z_{w}}{v_{w}\cos \theta }}-{\frac {2z_{w}\tan \theta }{v_{b}}}+{\frac {2x}{v_{b}}},}$ (C-49c)

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

 ${\displaystyle t_{-}=t_{+}+{\frac {2x}{v_{b}}}.}$ (C-49d)

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 ${\displaystyle t'_{ij}}$ for the refracted raypath from the shot location Sj to the receiver location Ri is given by

 ${\displaystyle t'_{ij}=t_{S_{j}B}+t_{BC}+t_{CR_{i}}.}$ (C-50)

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 θc = sin−1(vw/vb). 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.

By rewriting equation (C-50), for a flat or near-flat refractor, we obtain

 ${\displaystyle t'_{ij}={\frac {S_{j}B}{v_{w}}}+{\frac {DE-DB-CE}{v_{b}}}+{\frac {CR_{i}}{v_{w}}}.}$ (C-51)

By regrouping the terms, we get

 ${\displaystyle t'_{ij}=\left({\frac {S_{j}B}{v_{w}}}-{\frac {DB}{v_{b}}}\right)+\left({\frac {CR_{i}}{v_{w}}}-{\frac {CE}{v_{b}}}\right)+{\frac {DE}{v_{b}}}.}$ (C-52)

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

 ${\displaystyle t'_{ij}={\frac {z_{j}{\sqrt {v_{b}^{2}-v_{w}^{2}}}}{v_{b}v_{w}}}+{\frac {z_{i}{\sqrt {v_{b}^{2}-v_{w}^{2}}}}{v_{b}v_{w}}}+{\frac {x_{ij}}{v_{b}}}.}$ (C-53)

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:

 ${\displaystyle t'_{ij}=T_{j}+T_{i}+s_{b}x_{ij},}$ (C-54)

where

 ${\displaystyle T_{j}={\frac {z_{j}{\sqrt {v_{b}^{2}-v_{w}^{2}}}}{v_{b}v_{w}}},}$ (C-55)

 ${\displaystyle T_{i}={\frac {z_{i}{\sqrt {v_{b}^{2}-v_{w}^{2}}}}{v_{b}v_{w}}},}$ (C-56)

and

 ${\displaystyle s_{b}=1/v_{b}.}$ (C-57)

Note that Tj and Ti are the intercept time anomalies at shot and receiver locations, respectively, and sb is the bedrock slowness. Hence, for n shot/receiver stations the parameter vector is p : (T1, T2, …, Tn; sb). We refer to the scheme based on equation (C-54) as the variable-thickness GLI solution. Once the parameter vector p : (T1, T2, …, Tn; sb) is estimated, then the thickness of the weathering layer below shot and receiver locations can be computed using equations (C-55) and (C-56), respectively. For m picks of tij, and n + 1 parameters p : (T1, T2, …, Tn; sb), we have the following set of equations:

 ${\displaystyle {\begin{pmatrix}\vdots \\\vdots \\\vdots \\t'_{ij}\\\vdots \\\vdots \\\vdots \\\end{pmatrix}}={\begin{pmatrix}\\\\\cdots &1&\cdots &1&\cdots &x_{ij}\\\\\\\end{pmatrix}}{\begin{pmatrix}\vdots \\T_{j}\\\vdots \\T_{i}\\\vdots \\s_{b}\end{pmatrix}}.}$ (C-58)

By writing these equations in matrix notation, we have

 ${\displaystyle {\mathbf {t} '}={\mathbf {L} p},}$ (C-59)

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

 ${\displaystyle {\mathbf {e} }={\mathbf {t} }-{\mathbf {t} '}}$ (C-60)

is minimum and is given by:

 ${\displaystyle {\mathbf {p} }=({\mathbf {L} ^{T}L})^{-\mathbf {1} }{\mathbf {L} ^{T}t},}$ (C-61)

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:

1. Assume a value for the weathering velocity vw. This can be varied spatially based on available uphole information.
2. Estimate the parameter vector p : (T1, T2, …, Tn; sb), hence compute the intercept time anomalies at shot/receiver locations and the bedrock slowness by solving equation (C-61).
3. 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

 ${\displaystyle t'_{ij}=\alpha _{j}z_{j}+\alpha _{i}z_{i}+s_{b}x_{ij},}$ (C-62)

where

 ${\displaystyle \alpha _{j}={\sqrt {s_{wj}^{2}-s_{b}^{2}}},}$ (C-63)

 ${\displaystyle \alpha _{i}={\sqrt {s_{wi}^{2}-s_{b}^{2}}},}$ (C-64)

 ${\displaystyle s_{w}=1/v_{w}}$ (C-65)

and

 ${\displaystyle s_{b}=1/v_{b}.}$ (C-66)

Hence, for n shot/receiver stations, the parameter vector is p : (α1, α2, …, αn; sb). We refer to the scheme based on equation (C-62) as the variable-velocity GLI solution. Once the parameter vector p : (α1, α2, …, αn; sb) is estimated, then the weathering velocity below shot and receiver locations can be computed using equations (C-63) and (C-64), respectively. For m picks of tij, and n + 1 parameters p : (α1, α2, …, αn; sb), we have the following set of equations:

 ${\displaystyle {\begin{pmatrix}\vdots \\\vdots \\\vdots \\t'_{ij}\\\vdots \\\vdots \\\vdots \\\end{pmatrix}}={\begin{pmatrix}\\\\\cdots &z_{j}&\cdots &z_{i}&\cdots &x_{ij}\\\\\\\end{pmatrix}}{\begin{pmatrix}\vdots \\\alpha _{j}\\\vdots \\\alpha _{i}\\\vdots \\s_{b}\end{pmatrix}}.}$ (C-67)

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:

1. Specify a flat datum to which shot and receivers are to be lowered along vertical raypaths.
2. Estimate the parameters: p : (α1, α2, …, αn; sb) using the GLI solution given by equation (C-61) where the L matrix is defined as in equation (C-67).
3. 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 vw. 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]:

 ${\displaystyle t'_{ij}=s_{wj}Z_{j}+s_{wi}Z_{i}+s_{b}X_{ij},}$ (C-68)

where

 ${\displaystyle Z_{j}={\frac {z_{j}}{\cos \theta _{j}}},}$ (C-69a)

 ${\displaystyle Z_{i}={\frac {z_{i}}{\cos \theta _{i}}},}$ (C-69b)

 ${\displaystyle X_{ij}=x_{ij}-z_{j}\tan \theta _{j}-z_{i}\tan \theta _{i},}$ (C-69c)

 ${\displaystyle s_{w}={\frac {1}{v_{w}}},}$ (C-69d)

and

 ${\displaystyle s_{b}={\frac {1}{v_{b}}}.}$ (C-69e)

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

 ${\displaystyle {\begin{pmatrix}\vdots \\\vdots \\\vdots \\t'_{ij}\\\vdots \\\vdots \\\vdots \\\end{pmatrix}}={\begin{pmatrix}\\\\\cdots &z_{j}&\cdots &z_{i}&\cdots &X_{ij}\\\\\\\end{pmatrix}}{\begin{pmatrix}\vdots \\s_{wj}\\\vdots \\s_{wi}\\\vdots \\s_{b}\end{pmatrix}}.}$ (C-70)

Consider an initial estimate of the parameter vector p : (…, swj, …, swi, …, sb). 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 Δp in the paramater vector will change the modeled times as follows:

 ${\displaystyle \left[t'_{ij}\right]_{modeled}=\left[t'_{ij}\right]_{initial}+\left[{\frac {\partial t'_{ij}}{\partial p}}\right]_{modeled}\Delta p.}$ (C-71)

The error in modeling the traveltimes is given by

 ${\displaystyle e_{ij}=\left[t_{ij}\right]_{observed}-\left[t'_{ij}\right]_{modeled.}}$ (C-72)

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

 ${\displaystyle e_{ij}=\left[t_{ij}\right]_{observed}-\left[t'_{ij}\right]_{initial}-\left[{\frac {\partial t'_{ij}}{\partial p}}\right]_{modeled}\Delta p.}$ (C-73)

Now define the difference Δtij between the observed traveltimes tij and the initial estimate of the modeled traveltimes ${\displaystyle t'_{ij},}$ and rewrite equation (C-73) to get

 ${\displaystyle e_{ij}=\Delta t_{ij}-{\frac {\partial t'_{ij}}{\partial p}}\Delta p.}$ (C-74)

The second term on the right is the amount of change in Δtij as a result of the change in the parameter Δp. Define this term as ${\displaystyle \Delta t'_{ij},}$ and rewrite equation (C-74) once more to obtain

 ${\displaystyle e_{ij}=\Delta t_{ij}-\Delta t'_{ij},}$ (C-75)

where

 ${\displaystyle \Delta t'_{ij}={\frac {\partial t'_{ij}}{\partial p}}\Delta p.}$ (C-76)

The derivatives ${\displaystyle {\partial t'_{ij}/\partial p}}$ in equation (C-76) can be computed by differentiating the model equation (C-68) with respect to each of the parameters:

 ${\displaystyle {\frac {\partial t'_{ij}}{\partial s_{wj}}}\equiv Z_{j},}$ (C-77a)

 ${\displaystyle {\frac {\partial t'_{ij}}{\partial s_{wi}}}\equiv Z_{i},}$ (C-77b)

and

 ${\displaystyle {\frac {\partial t'_{ij}}{\partial s_{b}}}\equiv X_{ij}.}$ (C-77c)

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

 ${\displaystyle \Delta t'_{ij}=Z_{j}\Delta s_{wj}+Z_{i}\Delta s_{wi}+X_{ij}\Delta s_{b}.}$ (C-78)

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 : (Δsw1, Δsw2, …, Δswn, Δsb). We refer to the scheme based on equation (C-78) as the iterative GLI solution. For m picks of tij, and n + 1 parameters Δp : (Δsw1, Δsw2, …, Δswn, Δsb), we have the following set of equations:

 ${\displaystyle {\begin{pmatrix}\vdots \\\vdots \\\vdots \\\Delta t'_{ij}\\\vdots \\\vdots \\\vdots \\\end{pmatrix}}={\begin{pmatrix}\\\\\cdots &Z_{j}&\cdots &Z_{i}&\cdots &X_{ij}\\\\\\\end{pmatrix}}{\begin{pmatrix}\vdots \\\Delta s_{wj}\\\vdots \\\Delta s_{wi}\\\vdots \\\Delta s_{b}\end{pmatrix}}.}$ (C-79)

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

 ${\displaystyle \mathbf {\Delta t'} ={\frac {\partial {\mathbf {t} '}}{\partial {\mathbf {p} }}}\mathbf {\Delta p} }$ (C-80)

or

 ${\displaystyle \mathbf {\Delta t'} =\mathbf {L\Delta p} ,}$ (C-81)

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 ${\displaystyle {\partial t'/\partial p}}$ 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

 ${\displaystyle {\mathbf {e} }=\mathbf {\Delta t} -\mathbf {\Delta t'} }$ (C-82)

is minimum and is given by

 ${\displaystyle \mathbf {\Delta p} =({\mathbf {L} ^{T}L})^{\mathbf {-} 1}{\mathbf {L} ^{T}\ \Delta t},}$ (C-83)

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 : (sw1, sw2, …, swn, sb).
• Compute Δtij, the time difference between the picked (observed) times tij and the initial modeled times ${\displaystyle t'_{ij}}$.
• Estimate the change in parameters: Δp : (Δsw1, Δsw2, …, Δswn, Δsb), by way of the GLI solution given by equation (C-83).
• Update the paramater vector p + Δp, and compute new modeled times ${\displaystyle t'_{ij}}$.
• 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):

 ${\displaystyle t''_{ij}=t'_{ij}-s_{b}x_{ij},}$ (C-84)

so that

 ${\displaystyle t''_{ij}=T_{j}+T_{i}.}$ (C-85)

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:

 ${\displaystyle t''_{ij}=z_{j}{\sqrt {s_{wj}^{2}-s_{b}^{2}}}+z_{i}{\sqrt {s_{wi}^{2}-s_{b}^{2}}}.}$ (C-86)

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

 ${\displaystyle {\frac {\partial t''_{ij}}{\partial z_{j}}}={\sqrt {s_{wj}^{2}-s_{b}^{2}}}\equiv Z_{j},}$ (C-87a)

 ${\displaystyle {\frac {\partial t''_{ij}}{\partial z_{i}}}={\sqrt {s_{wi}^{2}-s_{b}^{2}}}\equiv Z_{i},}$ (C-87b)

 ${\displaystyle {\frac {\partial t''_{ij}}{\partial s_{wj}}}={\frac {z_{j}s_{wj}}{\sqrt {s_{wj}^{2}-s_{b}^{2}}}}\equiv S_{wj},}$ (C-87c)

and

 ${\displaystyle {\frac {\partial t''_{ij}}{\partial s_{wi}}}={\frac {z_{j}s_{wi}}{\sqrt {s_{wi}^{2}-s_{b}^{2}}}}\equiv S_{wi}.}$ (C-87d)

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

 ${\displaystyle \Delta t''={\frac {\partial t''}{\partial p}}\Delta p,}$ (C-88)

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

 ${\displaystyle \Delta t''=Z_{j}\Delta z_{j}+Z_{i}\Delta z_{i}+S_{wj}\Delta s_{wj}+S_{wi}\Delta s_{wi}.}$ (C-89)

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

 ${\displaystyle {\begin{pmatrix}\vdots \\\vdots \\\vdots \\\Delta t''_{ij}\\\vdots \\\vdots \\\vdots \\\end{pmatrix}}={\begin{pmatrix}\\\\\cdots &Z_{j}&\cdots &Z_{i}&\cdots &S_{wj}&\cdots &S_{wi}&\cdots \\\\\\\end{pmatrix}}{\begin{pmatrix}\vdots \\\Delta z_{j}\\\vdots \\\Delta z_{i}\\\vdots \\\Delta s_{wj}\\\vdots \\\Delta s_{wi}\\\vdots \\\end{pmatrix}}.}$ (C-90)

Follow the steps which involve equations (C-80) through (C-83) to estimate the parameter vector p : (…, zj, …, zi, …, swj, …, swi).

## C.10 L1-norm refraction statics

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

 ${\displaystyle E=\sum _{ij}\vert \vert t_{ij}-t'_{ij}\vert \vert ^{2},}$ (C-91a)

where tij are the actual traveltime picks and ${\displaystyle t'_{ij}}$ 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 L2 norm. For statics applications, it may be desirable to use the L1 minimization norm defined by

 ${\displaystyle E=\sum _{ij}\vert \vert t_{ij}-t'_{ij}\vert \vert .}$ (C-91b)

Outliers in picked times can cause biased results from the L2-norm, whereas, they may be better handled by the L1-norm schemes. The L2-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 L1-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 L1-norm refraction statics corrections.

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 L1-norm. Note that these results are comparable to those obtained from L2-norm minimization (Figure 3.4-17).

The results of the L1-norm statics estimates are summarized in Figure C-7. Compare with those of the L2-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 L1-norm solution to refraction statics is the sum of the differences between the observed picks tij and the modeled traveltimes ${\displaystyle t'_{ij}}$ (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 L1-norm solution. Again, results are comparable to those obtained from L2-norm minimization (Figure 3.4-26).

Results of the L1-norm refraction statics solution are summarized in Figure C-9. Compare with those of the L2-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 L1-norm solution is based on a single layer with a constant velocity of 1400 m/s as for the L2-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 L1-norm solution. The results are comparable to those obtained from L2-norm minimization (Figure 3.4-32).

## Equations

 ${\displaystyle t^{2}=C_{0}+C_{1}x^{2}+C_{2}x^{4}+C_{3}x^{6}+\cdots }$ (3)

 ${\displaystyle {\frac {\Delta f}{f}}={\frac {\Delta t_{NMO}}{t_{0}}}}$ (6)

 ${\displaystyle t^{2}=t_{0}^{2}+{\frac {x^{2}\sin ^{2}\alpha }{v^{2}}}}$ (7)

 ${\displaystyle t^{2}=t_{0}^{2}+{\frac {x^{2}\cos ^{2}\phi }{v^{2}}}}$ (10)

 ${\displaystyle v_{NMO}={\frac {v}{\cos \phi }}}$ (11)

 ${\displaystyle \varphi ={\frac {1}{2}}\left[\sin ^{-1}{\frac {v_{w}}{v_{b}^{-}}}-\sin ^{-1}{\frac {v_{w}}{v_{b}^{+}}}\right]}$ (46a)

 ${\displaystyle v_{b}={\frac {2\cos \varphi }{\begin{pmatrix}{\frac {1}{v_{b}^{-}}}+{\frac {1}{v_{b}^{+}}}\end{pmatrix}}}}$ (46b)

 ${\displaystyle z_{w}={\frac {v_{b}v_{w}t_{i}^{-}}{2\cos \varphi {\sqrt {v_{b}^{2}-v_{w}^{2}}}}}}$ (46c)

 ${\displaystyle t_{+}={\frac {2z_{w}{\sqrt {v_{b}^{2}-v_{w}^{2}}}}{v_{b}v_{w}}}}$ (48a)

 ${\displaystyle t_{-}=t_{+}+{\frac {2x}{v_{b}}},}$ (48c)

 ${\displaystyle T_{j}={\frac {z_{j}{\sqrt {v_{b}^{2}-v_{w}^{2}}}}{v_{b}v_{w}}}}$ (53a)

 ${\displaystyle T_{i}={\frac {z_{i}{\sqrt {v_{b}^{2}-v_{w}^{2}}}}{v_{b}v_{w}}}}$ (53b)

## References

1. 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. Castle (1994), Castle, R. J., 1994, A theory of normal moveout: Geophysics, 59, 983–999.
3. Slotnick, 1959, Slotnick, M. M., 1959, Lessons in seismic computing: Soc. Expl. Geophys.
4. Grant and West, 1965
5. De Bazelaire (1988), De Bazelaire, E., 1988, Normal moveout revisited: Inhomogeneous media and curved interfaces: Geophysics, 53, 143–157.
6. 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.
7. 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. 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. 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. 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. 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.
12. 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.
13. 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.
14. Farrell and Euwema, 1984, Farrell, R. C. and Euwema, R. N., 1984, Refraction Statics: Proc. Inst. Electr. Electron. Eng., 72, 1316–1329.
15. 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.