The Solution

With the creation of a number of functions for seismic file import, graphics display, and data analysis, Mathematica was converted into a powerful environment for well ties and wavelet extraction. The following sequence of inputs demonstrates an example workflow using this functionality in Mathematica 5.

The flow begins by loading the packages with the required utility. On production machines, this is typically autoloaded.

Next, we load in the data, consisting of the well log to be tied and a window of seismic data from receivers near the well location. In this case, we specify one set of data, or trace, from each side of the well location, forming a grid of traces. The data set can then be manipulated by a handle, which we return and here store in the variable dset.

Separating out the correlation stage from the data loading allows us to iterate correlations without the expense of reloading the data.

The well log data typically consists of a series of measurements of rock velocity and density at half-foot increments. Seismic data is measured in time, while the log data is in depth; these must be put in the same domain for further comparison. By finding the rock velocities within the well, we can estimate the travel time for a shock wave going from the surface, down to a particular depth, and back up to the receiver. This data is usually not complete, however, which may cause us to incorrectly estimate a travel time corresponding to depth. Well log data can also suffer from systematic errors in measurement, which can be accounted for when tying the log by applying a small time shift between the well log (sonic) and the seismic data. The user does this by specifying the correspondence between particular times in the two data sets. The sonic can also be stretched out in time to compensate for other systematic effects; this is given as a range of percentage stretches. A portion of the seismic and the sonic are then correlated over some range of lags, given below as 900 milliseconds of data correlated over ±200 milliseconds.

So, each extracted seismic trace is cross-correlated with each stretched sonic log. For the parameters used previously, this results in 18 time series of ±200 milliseconds in length. We are interested in knowing the value and location of the peak correlation coefficients. The function ShowCorrelations displays all the correlations, instantaneous amplitude envelopes, pointers to the peak locations, and so on. The bar on the left gives the color coding for the correlation pointers; this example underlays the data with a density plot of the amplitudes, which in some cases may make some features more obvious.

Examining the previous graph, we may decide that 1% stretch provides the best correlations. We select it for further operations.

The correlations are combined to give us a prototype for the wavelet. They may be combined in different ways, the simplest being a stack, or average. The ideal seismic wavelet in this context is zero-phase; for our purposes this means we want the phase spectra over the dominant frequency range to be near zero. ShowStretch can be used to display the characteristics of the correlations and the extracted wavelet (where the wavelet is the rightmost vertical trace in the following graphic).

Frequently, the integration processor will choose an appropriate linear phase transformation by hand to make the wavelet approximately zero-phase. However, automated techniques are also available. For example, here we perform a linear fit on the phase data, which can then be used to transform the phases via a linear phase shift.

We can view the transformed wavelet, which now looks much closer to the desired shape with a strong central peak and balanced side lobes.

This extracted wavelet and the well logs can then be used to generate a synthetic seismic trace, which can be compared to the actual seismic data to verify the quality of the tie. If the tie is good, we would expect events (significant peaks and troughs in the time series) to be aligned between the synthetic and the real seismic. The synthetic data is generated by assuming that a sufficiently large velocity or density contrast in the rock will correspond to a reflection; thus, from the well logs, we can generate a time series of reflection coefficients that, when convolved with a wavelet, should look much like the seismic data.

Seismic data is frequently collected from a single line of sensors, in which case the data is called 2D, or from a grid of sensors, where it is known as 3D. In this case, the data is 3D, but for conciseness we will only compare the synthetic to a single line of receivers via the Show3DFalse option. Furthermore, as indicated previously, the start of the well data is often questionable, so we frequently create an artificial well log for this starting section, which can be adjusted for the alignment we derived. Here we suggest to the routines that the first 800 feet of data can be replaced by an artificial curve.

Once we are satisfied with the results, the extracted wavelet, synthetic seismic, and modified logs can be saved for handling elsewhere.

Naturally, the previous flow demonstrates only a few of the many options available to a user; many of these functions include numerous options to control the nature of the display or the format of the data. A large selection of transformations to perform on the data are also present, and Mathematica's programming environment allows arbitrary additional manipulations as required.