Residual statics estimation by traveltime decomposition

Seismic Data Analysis

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

Static deviations from a hyperbolic traveltime trajectory are illustrated in Figure 3.3-17. After NMO correction, the misalignment of the wavelet associated with a reflection time horizon h along the offset axis in the CMP gather will yield a poor stack trace. We want to estimate the time shifts from the time of perfect alignment and correct for these time shifts. To do this, a model is needed for the moveout-corrected traveltime from a source location to a depth point on a reflecting horizon, then back to a receiver location. Figure 3.3-18 shows the geometry and notation that will be used in defining this model. The key assumption in the commonly used traveltime model discussed here is that residual statics are surface-consistent . This means that static shifts are time delays that depend solely on source or receiver locations at the surface, not on raypaths in the subsurface. Such an assumption is valid if all raypaths, regardless of source-receiver offset, are vertical within the near surface layer. The surface-consistent assumption usually is a good one, because the weathering layer often has very low velocity and a strong refraction at its base tends to make travel paths vertical. The assumption may not be good for a high-velocity permafrost layer that causes rays to bend away from vertical.

• 

  Figure 3.3-17  Picking traveltime deviations from NMO-corrected gathers.

• 

  Figure 3.3-18  Surface-consistent statics model to establish the traveltime model equation (25). Here, T = topographic layer, B = base of weathering layer, D = datum to which static corrections are made, R = deep reflector, j = shot station index, i = receiver station index, k = midpoint location index, x_ij = offset between the shot and receiver stations.

The picked traveltime t_ij that corresponds to the jth source station, the ith receiver station, and the kth midpoint [k = (i + j)/2] along a specific reflection time horizon h (Figure 3.3-17) can be approximately modeled by the following equation ; :

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

(25)

where s_j is the residual static time shift associated with the jth source station, r_i is the residual static time shift associated with the ith receiver station, and G_k is the difference in two-way time at a reference CMP location and the traveltime at the kth CMP location along the hth horizon. This term refers to structural variations along the horizon and is called the structure term. The term ${M_{k}}{x_{ij}^{2}}$ is the residual moveout which is assumed to be parabolic. It accounts for the imperfect moveout correction within the specific time gate that includes the reflection time horizon h. The coefficient M_k has the dimensions of (time / distance²).

To be specific in an analysis of the system of equations implied by equation (25), assume n_s shot locations, n_r receiver locations, and n_G CMP locations. Then define the fold as n_f. The purpose is to decompose the observed traveltimes estimated (picked) from the data t_ij to their individual components as defined on the right side of equation (25). The number of time picks (or individual equations) is equal to n_G × n_f. The number of unknowns is n_s + n_r + n_G + n_G. Typically, (n_G × n_f) > (n_s + n_r + n_G + n_G); hence, there are more equations than unknowns.

We can formulate a least-squares problem in which we must minimize the sum of the least-squares error energy between the observed traveltime picks t_ij and the modeled traveltimes $t'_{ij}$:

$E=\sum \limits _{ij}(t_{ij}-t'_{ij})^{2}.$

(26)

Residual statics corrections involve three phases:

1. Picking traveltime deviations t_ij based on crosscorrelation of traces in a CMP gather with a reference or pilot trace that needs to be defined in some fashion,
2. Modeling t_ij by way of equation (25) and decomposing it into its components: source and receiver statics, structural and residual moveout terms, and
3. Applying the derived source and receiver terms s_j and r_i, respectively, to traveltimes on the pre- NMO-corrected CMP gathers.

The picking phase relates to estimating traveltimes t_ij from the data. Several picking schemes are in use in the industry. The one described here often is known as a pilot trace scheme. Starting with the CMP gathers that are NMO-corrected using a preliminary velocity function(s), trace amplitudes are scaled to a common rms amplitude in the time gate(s) to be used for picking. For clarity, consider a time gate over the kth midpoint gather as illustrated in Figure 3.3-17. (It is preferable to start with a gather that has a good signal-to-noise ratio.) A stack trace then is constructed within the time gate h. Each individual trace in the gather is crosscorrelated with the stack trace. Time shifts t_ij, which correspond to maximum crosscorrelations, are picked. A preliminary pilot trace is constructed by stacking the time-shifted traces in the gather. This preliminary pilot trace is, in turn, crosscorrelated with the original traces in the gather and new values for t_ij are computed. A final pilot trace is constructed again by stacking the original traces shifted by the new values for t_ij. This final pilot trace is crosscorrelated with the traces of the next gather to construct the preliminary pilot trace for that gather. The process is performed this way on all CMP gathers moving to the left and right from the starting point. The picked time deviations t_ij are passed on to the next phase, which involves decomposing these time picks into components as defined by equation (25).

Several practical issues are involved in the picking phase. Band-pass filtering often helps to estimate a reliable time shift that corresponds to the peak crosscorrelation value. Selection of the time window used for crosscorrelation is another important factor. If needed, the window should be allowed to change laterally so that it follows the marker horizon(s). A maximum threshold for the correlation shift can be imposed to prevent unrealistically large time deviations from being passed on to the decomposition phase. Any time deviations greater than a specified maximum allowable shift can be set to that shift, or rejected altogether. Alternatively, the rejected shift can be replaced by a secondary correlation peak value. Finally, the input CMP gathers must be NMO-corrected by using a regional velocity function or velocities derived from a preliminary velocity analysis. These parameters are discussed in detail later in the section.

The next step in residual statics corrections involves least-squares decomposition of the time picks t_ij. For a basic understanding of this step, consider the following special problem. Suppose there are four observations t_ij measured at receiver locations x_i, where i = 1, 2, 3, 4. We want to fit the observed data into a straight line t = a + bx, which is best in a least-squares error sense. Start with the following set of equations:

${\begin{array}{lcr}t_{1}\approx a+bx_{1}\\t_{2}\approx a+bx_{2}\\t_{3}\approx a+bx_{3}\\t_{4}\approx a+bx_{4}.\\\end{array}}$

(27)

There are two unknowns, a and b, and four equations. This problem is similar to decomposing the time picks into various components as described by equation (25). We define the error series e_i, such that

${\begin{array}{lcr}-t_{1}+a+bx_{1}=e_{1}\\-t_{2}+a+bx_{1}=e_{2}\\-t_{3}+a+bx_{1}=e_{3}\\-t_{4}+a+bx_{1}=e_{4}.\\\end{array}}$

(28)

We want to minimize the energy of the cumulative errors as defined by

$E=\sum \limits _{i}e_{i}^{2}.$

(29a)

The energy for the ith error is computed by squaring both sides of equation (28)

$e_{i}^{2}=(-t_{i}+a+bx_{i})^{2}.$

(29b)

By substituting into equation (29a) and summing over i = 1, 2, 3, 4, we obtain

$E=\sum \limits _{i}t_{i}^{2}+4a^{2}+b^{2}\sum \limits _{i}x_{i}^{2}-2a\sum \limits _{i}t_{i}-2b\sum \limits _{i}x_{i}t_{i}+2ab\sum \limits _{i}x_{i}.$

(30)

To find the best line, we want to determine unknowns a and b so that error energy E is minimal. The sum E of the squared errors will attain a minimum if a and b are chosen so that

${\frac {\partial E}{\partial a}}=8a-2\sum \limits _{i}t_{i}+2b\sum \limits _{i}x_{i}=0,$

(31a)

and

${\frac {\partial E}{\partial b}}=2b\sum \limits _{i}x_{i}^{2}-2\sum \limits _{i}x_{i}t_{i}+2a\sum \limits _{i}x_{i}=0.$

(31b)

We now have two equations with two unknowns that can be put into a matrix form

${\begin{pmatrix}4&\sum \nolimits _{i}x_{i}\\\sum \nolimits _{i}x_{i}&\sum \nolimits _{i}x_{i}^{2}\end{pmatrix}}{\begin{pmatrix}a\\b\end{pmatrix}}={\begin{pmatrix}\sum \nolimits _{i}t_{i}\\\sum \nolimits _{i}x_{i}t_{i}\\\end{pmatrix}}.$

(32)

Equation (32) can be solved for the unknowns a and b. The minimum energy between the estimated model and the actual data values then can be computed by solving for a and b and substituting into equation (30).

Table 3-8. Time picks t_i at various x_i locations.

i	xi	$t'_{i}$
1	1	2.4
2	2	2.9
3	3	3.6
4	4	4.1

Consider the time values at various x locations specified in Table 3-8. Using these time picks, the elements of the coefficient matrix on the left side and the column matrix on the right side of equation (32) are

${\begin{pmatrix}4&10\\10&30\end{pmatrix}}{\begin{pmatrix}a\\b\end{pmatrix}}={\begin{pmatrix}13\\35.4\\\end{pmatrix}}.$

(33)

The solution for equation (33) is a = 1.8 and b = 0.58. Note that we are not constrained by the number of observations when setting up the least-squares problem. This problem involves more equations (27) than unknowns, a and b.

The least-squares approach has a wide range of applications in applied geophysics. We also formulate the inverse filtering used in deconvolution based on a least-squares procedure (Section B.5). The Wiener filter itself is based on a least-squares estimation of a future sample point in a given time series. Here, we discussed another application of the least-squares method — residual statics estimation by traveltime decomposition, which involves more equations than unknowns. In earth modeling in depth, when discussing 2-D surface data processing, we encounter the problem of fitting irregularly spaced observations into a uniform grid. This involves least-squares fitting to a local plane (Section J.5).

We now return to the problem of residual statics estimation and refer to equation (26). We have the observed traveltime deviations t_ij picked from moveout-corrected CMP gathers within a specified time gate that includes a strong reflection. We want to estimate the parameters — the surface-consistent source and receiver static shifts, s_j and r_i, respectively. These parameters are related to the modeled traveltime deviations $t'_{ij}$ by way of the model equation (25). If there are 50 000 picks, then you have 50 000 model equations.

Substitute for t_ij from equation (25) and minimize the error energy E in equation (26) by requiring

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

(34)

which yields n_s + n_r + n_G + n_G equations and that many unknowns.

These equations can be solved for the residual statics associated with n_s source locations, n_r receiver locations, n_G structural terms, and n_G residual moveout terms. For a real data set, the number of these terms can be large. It is easy to solve the problem with the two parameters given by equation (34). However, when dealing with a large number of linear equations, the solution must be done accurately and efficiently. Wiggins et al. (1978) used the Gauss-Seidel iterative procedure to solve for equations (34) (Section C.4).

The Gauss-Seidel method is best described by returning to the earlier line fitting example and solving for a and b in equation (33). When written as normal equations, we have

$4a+10b=13$

(35a)

and

$10a+30b=35.4.$

(35b)

These equations are rearranged as

$a=3.25-2.5b$

(36a)

and

$b=1.18-0.333a.$

(36b)

Since the Gauss-Seidel technique is iterative, starting values are needed. To start the iteration, set a = b = 0. Substituting b = 0 in equation (36a), we obtain a = 3.25. Putting this updated solution into equation (36b) gives b = 0.0977. At the end of the first iteration, a = 3.25 and b = 0.0977. For the second iteration, put b = 0.0977 into equation (36a) to get a = 3.0075. Put a = 3.0075 into equation (36b) to get b = 0.1785. This iterative procedure is continued in Table 3-9.

The solution from the iterative procedure slowly converges toward the actual values, a = 1.8, b = 0.58. Convergence is not always guaranteed. Nevertheless, convergence can be attained with the Gauss-Seidel method, provided the unknowns are ordered properly; that is, iteration starts with the correct unknown. This problem is addressed in Exercises 3-11 and 3-12. The advantage of the Gauss-Seidel method is its ability to solve the large number of simultaneous equations rapidly.

Table 3-9. Gauss-Seidel iteration for solving equation (33) for a and b. (Values were rounded off for tabulation.)

Actual values: a = 1.8, b = 0.58.
Iteration	a	b
1	3.25	0.0977
5	2.4918	0.3502
10	2.0712	0.4902
15	1.9031	0.5462
20	1.8358	0.5686
25	1.8090	0.5779
29	1.7997	0.5807

The question of when to terminate the iterations remains. The rate at which the solution changes can be examined after each iteration. The computation is stopped when this rate falls below a specified threshold value.

The starting values for solving the normal equations (34), can be chosen as s_j = r_i = G_k = M_k = 0, as in the simple numerical example. To some extent, the best order of iteration depends on the exact nature of the statics problem at hand. One order of iteration that is commonly in use follows: Compute the structure term G, the residual moveout term M, the static shift associated with source location s, then the static shift associated with receiver location r. The procedure cycles back to G in the next iteration and continues until convergence is satisfactory. The order in which the individual terms are computed theoretically can be interchanged. However, the above order forces long-wavelength variations of the time picks to be concentrated mainly in the structure term. This leads to a lesser number of iterations (typically 2 or 3) for wavelength components of the statics that are less than half the largest data offset. A large number of iterations are required to handle static variations that are greater than the maximum offset.

After computing the individual static shifts associated with each source and receiver location, these shifts are passed on to the application phase, whereby the shifts are applied to the pre-NMO-corrected gather traces (Figure 3.3-13). In areas with extremely poor signal-to-noise ratio or complicated near-surface variations, multiple passes of residual statics corrections may be necessary. In other words, output from the first pass may be NMO-corrected, new picks estimated, decomposed, and applied, and so on.

Figure 3.3-13  Moveout-corrected common-shot gathers from the same land line as in Figure 3.3-1. Common-shot gathers (a) before and (b) after residual statics corrections.

The question remains as to whether equations (25) are independent of each other. This is important, for it tells whether the solution is unique or not. For simplicity, consider the case of zero structure and zero residual moveout. In that case, equations (25) take the form

$t_{ij}=s_{j}+r_{i}.$

(37a)

Consider only four unknowns: s_j, where j = 1, 2, and r_i, where i = 1, 2, which yield four equations:

${\begin{array}{lcr}t_{11}=s_{1}+r_{1}\\t_{12}=s_{2}+r_{1}\\t_{21}=s_{1}+r_{2}\\t_{22}=s_{2}+r_{2}.\\\end{array}}$

(37b)

However, from a close examination of equations (37b), we note that

$t_{11}+t_{22}=t_{12}+t_{21}.$

(38)

Therefore, one of the four equations (37b) is redundant, leaving three independent equations for four unknowns.

This simple exercise shows that the statics solution obtained by solving equations (34) suffers from an uncertainty. In particular, the solution obtained by, say, the Gauss-Seidel iteration, is nonunique. It is but one of many possible solutions. In fact, the solution may not even be physically reasonable. Because of the problem of fewer independent equations than unknowns, it may be necessary to impose a constraint on the solution from the traveltime decomposition (equation 25):

1. A plausable constraint is that the difference between shot and receiver statics be minimal . It can be argued that this may not be valid, as is the case with dynamite data in which shots and receivers do not occupy the same physical locations.
2. Other possible constraints can be in the form of restrictions on the spatial variation of structure, moveout, or statics terms themselves; all are used in various practical implementations. For instance, one may postulate a model equation for the residual statics problem, where the structure term has been omitted. Then, you have to design your picking strategy accordingly, such that you do not accumulate the traveltimes as you go from one CMP location to another. This may be preferable in areas with very low-relief structures.
3. We may want to impose the surface-consistency rule in a strict sense, and assign the same static shift to a shot and receiver if they occupy nearly the same physical location.
4. We may also opt to set the residual moveout coefficient M to a constant across the line if we are fairly confident of the velocity control. This then gives us fewer parameters to estimate, thus making the residual statics solution more likely to be stable and physically plausable.

See also

• Residual statics estimation by stack-power maximization
• Traveltime decomposition in practice
• Maximum allowable shift
• Correlation window
• Other considerations
• Stack-power maximization in practice
• Exercises
• Topics in moveout and statics corrections

External links

find literature about Residual statics estimation by traveltime decomposition