Gaussian-process depth correction

1Priorbefore any well is seen

Before it has seen a single well, the Gaussian process places a probability distribution over what the correction surface might look like. It needs only one rule: the covariance of r at any two points depends only on the distance d between them, and is given by the kernel:

k(d) = σ_f² exp(−d² / 2ℓ²) + σ_n² [d = 0]

Each of the three hyperparameters controls one thing. The signal standard deviation σ_f sets the amplitude of the correction surface. The length scale ℓ sets how quickly the correlation decays with distance: at a separation of ℓ the correlation coefficient is e^−1/2 ≈ 0.61, and at 2ℓ only 0.14 remains. The noise standard deviation σ_n appears only at d = 0: it is the part that is unique to each well and cannot be predicted from its neighbours; tie errors, the pick of the log top and local velocity anomalies all end up here.

With normalize_y=True, sklearn first standardises r: it subtracts the mean at the wells, 21.84 m, and divides by the population standard deviation, 153.41 m. The fitted hyperparameters are all on this standardised scale: amp2 = 1.188, ℓ = 5.58 km, noise = 0.0883; converted back to metres, σ_f = √amp2 × 153.41 = 167 m and σ_n = √noise × 153.41 = 45.6 m. The prior mean is the mean at the wells, so far from all wells the correction returns to the constant 21.84 m. When some wells are switched off, sklearn re-standardises with the mean and standard deviation of the remaining wells, while the hyperparameters stay in standardised units; the metre values on the sliders are therefore always converted with the 153.41 m of the full well set, and the actual values in metres change slightly with the set of active wells.

2Conditioningon the calibration wells

Bring in the wells and the prior becomes the posterior: of all possible surfaces, only those that pass close to the wells, with misfits in proportion to the noise, remain. For a Gaussian process this step has a closed form. Let ȳ be the mean of r over the active wells and s its population standard deviation, z = (y − ȳ) / s the standardised observations, K the well-to-well kernel matrix and k_* the kernel values between the target point and each well (all in standardised units):

r̂(x_*) = ȳ + s · k_*^T (K + noise · I)⁻¹ z

sd_latent = s · √( amp2 − k_*^T (K + noise · I)⁻¹ k_* )

sd_pred = s · √( amp2 + noise − k_*^T (K + noise · I)⁻¹ k_* )

The page solves (K + noise · I) by Cholesky decomposition and, like sklearn, adds a tiny alpha to the diagonal. The latent sd is the uncertainty of the correction surface itself; the predictive sd adds the white noise and describes where r at a new well might fall. For this kernel, sklearn's predict(return_std=True) returns the latter.

The posterior surface is two-dimensional; Fig. 2 is only a slice through it along one line. Wells seldom lie exactly on the line, so the figure shows the wells within 1.5 km of it, projected onto the line (the farther away, the paler), and the curve need not pass through them: the wells lie off the line, and the noise term allows the curve to pass beside a well in any case. There are currently 52 active wells, standardised with ȳ = 21.84 m and s = 153.41 m.

3Whole surveythe correction map

The same calculation over the whole survey: one point every 16 traces (spacing 0.4 km), 166 × 204 grid points in all, drawn only where the supplied surface is defined. The left map is the correction r̂: blue where the V(z) depth is too shallow (r̂ < 0, deepened by the correction), red where it is too deep (made shallower). The right map is the posterior standard deviation: smallest within the well clusters, it approaches the prior value away from the wells, which is 173 m for the predictive sd and 167 m for the latent sd.

Coordinates x and y are distances in km from a corner of the seismic survey. On the maps x increases to the right and y downwards, so the top of each map faces south-west: y increases towards an azimuth of 49° and x towards 319° (the arrow at top right points north).

4Leave-one-well-outeach well predicted by the others

The misfit at the wells says little: the posterior mean is drawn towards the wells by construction. The real test is leave-one-well-out: remove one well, redo normalize_y and the conditioning on the remaining 51 with the same hyperparameters, and predict the well that was removed, exactly as in the original sklearn analysis. Wells that are switched off take no part in the fit; their errors are computed directly from the posterior mean of the active wells, so they are already held out.

With the fitted parameters, the median absolute error at the 52 wells falls from 107 m uncorrected to 37 m, and the RMS from 155 m to 112 m. Drag the sliders and these numbers are recomputed.

The right panel plots the median error as a function of the length scale (σ_f and σ_n at their current values). When ℓ is very short, almost no two wells are correlated and the well left out receives only the mean of the others, which amounts to a constant correction; when ℓ is very long, the whole surface is smeared into a gentle slope. The fitted value lies in the trough between the two. The floor of the trough is not smooth: the median is set by the one or two wells in the middle of the ranking and jumps by a few metres when ℓ changes slightly; and the marginal likelihood measures the overall probability, not the median of these 52 numbers, so the lowest point of the trough need not coincide with the fitted value. Of the 24 sample points on the curve, the lowest is at ℓ = 7.55 km (37.5 m); at the current ℓ (orange dot) the median is 36.6 m.

5Marginal likelihoodwhere the fitted parameters come from

The fitted values come from maximising the log marginal likelihood. For the standardised observations z (n wells):

log p(z | θ) = −½ z^T (K + noise · I)⁻¹ z − ½ log |K + noise · I| − (n/2) log 2π

The first term is the data fit: the better the data at the wells match the correlation structure that the model assumes, the larger it is. The second term is a complexity penalty: the shorter ℓ, the less correlated the wells, the more surfaces the model can generate, the larger the determinant |K + noise · I| and the more negative this term; conversely, the longer ℓ and the smaller the noise, the more rigid the model and the smaller the determinant, but as soon as the data at the wells are less smooth than the model assumes, the first term quickly becomes negative. The maximum is a compromise between the two. sklearn maximises it with L-BFGS-B within the bounds (ℓ ∈ [0.3, 60] km, amp2 ∈ [0.01, 100], noise ∈ [0.001, 10]), starting from the initial values and several random starting points, and reaches −54.34.

Fig. 5 is a slice through the ℓ × σ_n plane, with σ_f fixed at the current slider value (by default the fitted value); the colour shows how many log units below the maximum within the slice each point lies. The log marginal likelihood at the current parameters is −54.34. Switching wells off changes this surface, and the original fit is then no longer the optimum; ‘Re-maximise’ repeats the fit on the currently active wells (the page uses multi-start Nelder–Mead with the same bounds as sklearn).

6What it cannot fixlimits of a velocity correction

This correction fixes velocity only. It assumes that the supplied time surface is correct at the wells, which is why r was estimated only from wells whose tie agrees with the surface to within 50 ms. Another 32 wells took no part in the fit; they fall into two groups:

• 12 independent check wells: their tie is level C, but the surface agrees with it, so their r still reflects velocity and they provide a genuinely independent test. The correction reduces the median absolute error from 121 m to 43 m and raises the proportion of errors within ±100 m from 25 % to 83 %.
• 20 wells where the surface and the tie differ by more than 50 ms: at least one of the two is wrong, and many level-C ties are themselves derived from the regional V(z). Here r mainly reflects a difference in the time of the surface rather than velocity, and the correction only reduces the median absolute error from 492 m to 394 m.

The reason is simple: r̂ is a smooth surface, and where the surface is defined its 5th–95th percentile range spans only −85 m to 266 m, whereas the errors at these wells reach several hundred metres and change abruptly: at 48/11b-10 and 48/11b-6, not far apart, they are −738 m and 585 m. This is a problem of the surface, beyond the reach of a velocity correction. Take 48/02b-3 (a level-B tie): the supplied surface is 308.6 ms later than the tied top, so the converted depth is 752.0 m deeper than the log top (785.5 m TVDSS), while the correction there is only 85.1 m, leaving 666.9 m.