# The least-squares method

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

We want to estimate the near-surface parameters — weathering and bedrock velocities and thickness of the weathering layer at shot-receiver locations by least-squares inversion of the observed (picked) refracted arrivals. Formulation of this problem using the least-squares inversion leads to an estimate of the nearsurface parameters such that the difference between the observed arrivals and the modeled refracted arrivals is minimum in the least-squares sense. This method is not only applicable to 2-D line shooting but also to 3-D swath shooting geometries.

There are several ways to parameterize the nearsurface layer. The most general formulation would include varying weathering and bedrock velocities and the thickness of the weathering layer at all shot-receiver stations. This, however, would require linearizing the problem and iterating over the estimated parameters. The problem also would have to be constrained to stabilize the inversion. In a simplified version of this general formulation, weathering velocity may be fixed and assumed to be known. This leaves the weathering thickness and bedrock velocity as spatially varying parameters.

As for any inversion problem, we shall need a model equation that relates the model parameters we want to estimate to the modeled refracted arrival times. Refer to the sketch of a near-surface model in Figure 3.4-13. If we let the weathering thickness vary as might be the case in many field data applications, then we have a problem of not being able to write down an analytic expression for the refracted raypath. We would not even begin to anticipate how the head wave would behave or develop in case of a laterally varying refractor geometry.

Instead, to be able to make use of the first-break picks, we shall take a simpler approach. We want to describe the near-surface with minimal parameterization and consider the model with a flat refractor. Now, we can express the modeled traveltime ${\displaystyle t'_{ij}}$ for the critically refracted raypath from the source location Sj to the receiver location Ri (Figure 3.4-13) as

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

Keep in mind that we had to make the same assumption — that the refractor is flat or nearly-flat within the spread length, as in the generalized reciprocal method. The first and the third terms in equation (51a) are associated with the raypaths within the weathering layer, and the second term is associated with the raypath within the bedrock along the refractor. When the refractor dip is taken into account, the problem cannot be readily linearized.

By regrouping the terms in equation (51a), we get the expression

 ${\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}}}.}$ (51b)

Finally, rewriting in terms of the model parameters that we want to estimate — vw, vb, and zw, we obtain the expression for 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}}}.}$ (51c)

In addition to assuming a flat refractor, we fix the bedrock velocity vb but retain it as a parameter to be estimated. We also assume that vw is known. Under these assumptions, the model equation (51c) for the refracted arrivals can be written in the form [1]

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

where sb = 1/vb is the bedrock slowness, and

 ${\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)

By comparing equations (53a,48b), note that Tj and Ti actually are (half) intercept time values at the shot and receiver locations. Hence, for n shot-receiver stations, the parameter vector contains (T1, T2, …, Tn; sb).

 ${\displaystyle t_{i}={\frac {2z_{w}{\sqrt {v_{b}^{2}-v_{w}^{2}}}}{v_{b}v_{w}}},}$ (48b)
Figure 3.4-13  Geometry of refracted arrival used in deriving the least-squares solution for intercept times. Here, Sj and Ri are source and receiver stations, respectively; θc is the critical angle of refraction, zj and zi are depths to the bedrock at source and receiver locations, and vw and vb are weathering and bedrock velocities, respectively.

Equation (52), which models refracted arrivals, is solved for the parameter vector in the same manner using the generalized linear inversion (GLI) theory (Section C.8) as for equation (25), which models traveltime deviations associated with residual statics (Section C.4). The generalized linear inversion solution is based on the objective of minimizing the least-squares difference between the observed refracted arrival times tij and the modeled times ${\displaystyle t'_{ij}}$ defined by equation (52). Alternatively, the refraction statics solution can be obtained by using the minimization criterion that is based on the L1 norm (Section C.10).

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

An extension of equation (52) to represent the refraction statics associated with each shot and receiver location in two parts — long- and short-wavelength components, is proposed by [2]. This is done by including two more terms in equation (52) to represent the shot and receiver short-wavelength variations in refraction statics caused by rapid changes in elevation and near-surface layer geometry. The conjecture for such an extension is to stabilize the long-wavelength solution derived from equation (52). Albeit the extension includes short-wavelength terms for shots and receivers, you still would need to estimate residual statics and apply them to your data (equation 25).

Equation (52) describes a variable-thickness scheme since the weathering velocity is assumed to be constant. Note that the traveltime model equation (52) suffers from uncertainty in the value for the weathering layer velocity, as was the case for the generalized reciprocal method. The estimated refractor shape using the variable-thickness scheme (whether it is based on the the generalized reciprocal or least-squares method) does not yield the true refractor shape. Instead, the uncertainty in weathering velocity significantly influences the implied refractor shape. Perhaps, uphole information can be used to calibrate the estimated thicknesses and thus an acceptable weathering velocity can be chosen accordingly.

Once the parameter vector is estimated, then the thickness of the weathering layer below shot and receiver locations can be computed using the expressions for the intercept time values (equations 53a, 53b).

Because equation (52) does not contain a structure term, any long-wavelength static anomaly is partitioned between the other terms. This is not the case for the reflection-based residual statics model that is based on equation (25). Thus, following the field statics corrections (to account for elevation changes), statics corrections are estimated and applied in two stages:

1. Refraction-based statics corrections to remove long-wavelength anomalies, and
2. Reflection-based residual statics corrections to remove any remaining short-wavelength static shifts.

Both the generalized reciprocal and least-squares methods are based on computing intercept time anomalies at shot-receiver stations. The generalized reciprocal method yields multi-valued intercept times for each shot-receiver station, which need to be reduced to a single-valued intercept time profile along the line before computing the shot-receiver statics. On the other hand, the least-squares method yields a unique intercept time value for each shot-receiver station based on the least-squares minimization.

The generalized reciprocal method requires a special combination of raypaths associated with the traveltime picks that correspond to the terms in equations (50a, 50b) in estimating intercept times. These raypaths are not always attainable by 3-D swath shooting, and therefore, the generalized reciprocal method is most suitable for 2-D seismic profiles. On the other hand, the least-squares method yields intercept times for arbitrary shot-receiver locations associated with 2-D or 3-D recording geometries. The important point to keep in mind, however, is that neither method yields the true physical parameters for the near-surface — the variablethickness solution depends on the assumed value for the weathering velocity. However, with extra information such as from upholes, these solutions may be calibrated.