Cross-balancing intersecting surveys

William S. Harlan

September, 2003

Introduction

Tying seismic sections with well picks typically requires estimating a source wavelet from a sonic log. Few want to repeat this effort for each dataset with different recording parameters. Each frequency of a source wavelet can have a different amplitude and phase. A typical phase rotation or amplitude scaling does not change with frequency and cannot make one arbitrary wavelet look like another. Instead we should estimate different wavelets directly from cross-correlations of all intersecting datasets.

A collection of independent seismic surveys, 2D lines or 4D vintages, intersect each other, but reflections do not overlay well. Seismic attributes differ and complicate any stratigraphic interpretation, particularly of thin beds. We want to apply wavelet corrections to each of these datasets so that intersections tie as well as possible.

One approach was independently implemented by David Mackidd of EnCana, and by Bishop and Nunns [1]. This approach cross-correlates traces from different surveys near intersection points. From these cross-correlations are extracted time shifts, phase rotations, and scale factors by least-squares means. Bulk phase-rotations require care with phase-unwrapping, and time shifts are not independent of phase. These limited corrections may be robust in the presence of defective spectra, but at the cost of ignoring arbitrary phase and amplitude changes with frequency.

Instead, we take the approach of Henry and Mellman [3], and invert cross-correlations directly for arbitrary wavelets. This problem is very similar to the more-familiar problem of surface-consistent deconvolution (Levin [4]), which solves for prediction-error filters. I have implemented a variation on both. I solve the same wavelet equalization problem as Henry and Mellman, but with a non-quadratic objective function more like surface-consistent deconvolution, and with the Gauss-Newton optimization of Harlan [2].

Cross-correlations

Let us index independent seismic lines or datasets by indices $i$ and $j$. Certain pairs of lines cross at intersections indexed by $m$. Each $m$ will index a triplet of indices $\{m,i,j\}$. At each of these intersections $m$, a representative trace $d_i^m(t)$ from line $i$ will be cross-correlated with a trace $d_j^m(t)$ from line $j$, all as a function of time $t$. Let us define the cross-correlation $c_{ij}^m(t)$ as a function of lag time $t$ by

$$c_{ij}^m(t) = \int d_i^m(\tau-t) d_j^m(\tau) d \tau$$ (1)

$$= \int d_i^m(\tau) d_j^m(\tau+t) d \tau$$ (2)

$$= d_i^m(t) \star d_j^m(t)$$ (3)

$$= d_i^m(-t) \ast d_j^m(t) .$$ (4)

Here the star $\star$ is a concise notation for correlation, and the asterisk $\ast$ for convolution.

Typically cross-correlations will be averaged from multiple traces, but I will assume a single trace in the derivation of a solution.

The model

Assume that a data trace $d_i^m(t)$ is a convolution of a short wavelet $w_i (t)$ that is specific to the line $i$ with a reflectivity $r^m (t)$ that is specific to the intersection point $m$. Also assume additive noise $n_i^m (t)$ that is specific to the trace:

$$d_i^m(t) = \int r^m(t-\tau) w_i(\tau) d \tau + n_i^m(t)$$ (5)

$$= r^m(t) \ast w_i(t) + n_i^m(t) .$$ (6)

Assume the reflectivity is white:

$$r^m(t) \star r^m(t) = R^m \delta(t) ,$$ (7)

where $R^m$ is a constant, and $\delta (t)$ is an impulse function.

Assume noise $n_i^m (t)$ is uncorrelated with anything else:

$$n_i^m(t) \star n_i^m(t) \equiv N (t).$$ (8)

$$n_i^m(t) \star n_j^m(t) = 0, \ \forall i \neq j .$$ (9)

$$n_i^m(t) \star r^m(t) = 0, \ \forall i.$$ (10)

$$n_i^m(t) \star w_j(t) = 0, \ \forall i,j .$$ (11)

A correlation of two traces can then be rewritten as

$$c_{ij}^m(t) = d_i^m(t) \star d_j^m(t)$$ (12)

$$= [r^m(t) \ast w_i(t) + n_i^m(t)] \star [r^m(t) \ast w_j(t) + n_j^m(t)]$$ (13)

$$= [r^m(-t) \ast w_i(-t) + n_i^m(-t)] \ast [r^m(t) \ast w_j(t) + n_j^m(t)]$$ (14)

$$= [r^m(-t) \ast w_i(-t) \ast r^m(t) \ast w_j(t) ] + [n_i^m(-t) \ast n_j^m(t)]$$ (15)

$$= [r^m(t) \star r^m(t)] \ast [w_i(t)\star w_j(t) ] + [n_i^m(t) \star n_j^m(t)]$$ (16)

$$= R^m [ w_i(t) \star w_j(t) ] + \delta_{i-j} N (t).$$ (17)

We will only examine cross correlations where $i \ne j$, so the discrete delta function $\delta_{i-j}$ will be 0.

The objective function

Thus far the assumptions are very similar to those of Henry and Mellman [3]. They use time shifts as separate parameters so that their short wavelets are centered around zero-lag. Their objective function also measures differences between third-order convolutions of measured correlations with wavelets. These higher-order terms are easier to optimize with least-squares, but distribute errors more unpredictably.

I prefer to use longer wavelets to avoid a separate parameterization of time shifts. I directly minimize errors in cross-correlations modeled from the estimated wavelets. The objective function finds the collection of wavelets $w_i (t)$ and scale factors $R^m$ that best minimize

$$\min_ {w_i(t), R^m } \sum_{\{m,i,j\}} \int \{ c_{ij}^m(t) - R^m [ w_i(t) \star w_j(t) ] \}^2 dt + \epsilon \sum_i \int [w_i (t)]^2 dt .$$ (18)

Sum over the triplets of indices. The damping suppresses unnecessary frequencies that do not contribute significantly to the cross-correlations. The damping factor $\epsilon$ is an appropriate ratio of expected variances for noise to that of wavelets. Conservatively small values are sufficient for stable inverses. I set the variance of a wavelet sample to 1000 times the variance of a sample of noise, but this ratio can vary orders of magnitude without significantly affecting the result.

If reflectivities have reasonably consistent strengths, then one can safely set all scale factors $R^m$ to 1. This implementation makes this assumption.

Optimization

This objective function is a non-quadratic function of the wavelets $w_i (t)$, but it is very amenable to a Gauss-Newton algorithm that iteratively approximates the objective function as a quadratic.

Initialize all wavelets to delta functions. Linearize a perturbation of modeled cross-correlations $\Delta c_{ij}^m(t)$ with respect to perturbed wavelets $\Delta w_i(t)$:

$$\Delta c_{ij}^m(t) \equiv R^m [ w_i(t) + \Delta w_i(t)] \star [w_j(t) + \Delta w_j(t) ]$$ (19)

$$- R^m [ w_i(t) \star w_j(t) ]$$ (20)

$$\approx R^m [ w_i(t) \star \Delta w_j(t) + \Delta w_i(t) \star w_j(t) ] .$$ (21)

With this linearization, and its linear adjoint, the objective function becomes a least-squares (quadratic) function of the wavelet perturbations. We can use a standard conjugate-gradient algorithm to solve for the wavelet perturbations that best fit the correlations not yet modeled by reference wavelets. Perturbations are scaled appropriately (by performing a line search on the original objective function) before adding to the reference wavelets.

Correction of data with estimated wavelets

We can imagine using estimated wavelets directly for least-squares deconvolution of each line. This would however unnecessarily attempt to whiten the frequency spectra. Most likely, the data are already whitened as much as desired. We also do not want to compute inverse wavelets explicitly, which would amplify some otherwise very weak noisy frequencies.

To avoid unnecessary alteration of the spectral color, we want to reconvolve all deconvolved traces by a wavelet for a preferred reference line. This reference line may also be one that has been carefully tied to sonic logs or formation tops by synthetic seismograms. If more than one line is constrained, then I share the same wavelet for all during optimization. The estimated wavelet for these constrained lines is set aside as the reference wavelet.

For each of the unconstrained lines I solve for a transfer function mapping the estimated wavelet to the reference wavelet, again with least-squares damping. Finally the estimated transfer functions are convolved with the original data, combining the effects of deconvolution and reconvolution.

We can also use the estimated wavelets for each line to extract simpler corrections, such as a bulk time shift, instantaneous phase rotation, and a scale factor, as did Bishop and Nunns [1]. These limited corrections should be more stable and less sensitive to noise. Extracting phase corrections from estimated wavelets is also more robust than fitting phase shifts picked from the data. But for any frequencies shared by intersecting lines, we know we can do better and detect more arbitrary differences.

It is possible that certain lines will contain frequencies not present in any intersecting lines. These isolated frequencies will be absent from cross-correlations. A damped least-squares fit to cross-correlations will omit these frequencies from the estimated wavelet. To retain these frequencies, we can include autocorrelations as well as cross-correlations in the optimizations. A final correction to a reference wavelet may remove the isolated frequencies again.

References

1

Thomas N. Bishop and Alan G. Nunns.
Correcting amplitude, time, and phase mis–ties in seismic data.
Geophysics, 59:946–953, 1994.

2

William S. Harlan.
Simultaneous velocity filtering of hyperbolic reflections and balancing of offset–dependent wavelets.
Geophysics, 54(11):1455–1465, 1989.

3

M. Henry and G. R. Mellman.
Linearized simultaneous inversion for source wavelet equalization and mis–tie adjustment.
In 58th Annual International Meeting, SEG, Expanded Abstracts, pages 953–955. Soc. of Expl. Geophys., 1988.

4

S.A. Levin.
Surface–consistent deconvolution.
Geophysics, 54(09):1123–1133, 1989.