Correlation Length and Fractal Dimension Interpretation from Seismic Data Using Variograms and Power Spectra


Ken Mela and John Louie

Nevada Seismological Laboratory 174, University of Nevada, Reno, NV 89557

goldprosp@yahoo.com; louie@seismo.unr.edu

Click here for a 1.7 Mb PDF version including figures.


ABSTRACT

Recent modeling techniques for fluid flow in a reservoir or an aquifer require characterizing the statistics of the spatial variability of physical properties such as porosity or reflectivity. The major interest of such modeling is in horizontal variations. We acquired a non-traditional seismic survey and analyzed it for the statistical parameters of correlation length and fractal dimension. Variograms and power spectra can extract these parameters from migrated sections. Statistical parameters are valuable for newly established hydrogeologic modeling techniques, with their definition greatly increasing the understanding of the aquifer under study. Applying similar techniques may also prove valuable for reservoir exploitation in the petroleum industry. We acquired our survey on top of a bench of an open-pit diatomite quarry and extracted the statistical parameters for an exposed vertical face. The imaged zone of interest is therefore physically accessible to future permeability studies that could validate our interpreted parameters. We can interpret horizontal, but not vertical, correlation lengths and fractal dimensions from either variograms or power spectra of our migrated seismic images. The horizontal correlation length and fractal dimension extracted from the seismic variograms match those extracted from photographs taken of the quarry face.



INTRODUCTION

In July 1996 we acquired four high-resolution seismic lines on a nearly level bench of Eagle Picher's section 8 diatomite quarry south of Hazen, Nevada (approximately 40 miles east of Reno). We processed the data from these lines and calculated variograms and power spectra in horizontal and vertical directions from migrated seismic reflectivity sections. A classical variogram, and a moving window variogram estimator after Li and Lake (1994), were calculated for each section. From these processes we interpret horizontal correlation length and fractal dimension for the Tertiary lithology sampled.

Correlation length is the distance from a point beyond which there is no further correlation of a physical property associated with that point. Values for a given property at distances beyond the correlation length can be considered purely random. The horizontal correlation length is needed for determining the macrodispersivity (a parameter that indicates the distribution of concentrations of a contaminant front) used in stochastic contaminant transport and fluid flow models (Gelhar et al., 1979).

Fractal dimension is a parameter that can be used to determine the tortuosity of flow; i.e., the greater the fractal dimension, the farther a fluid volume must flow to travel between two points, or the less linear the flow between these two points will be. The fractal dimension can be used in simulations for design of fluid extraction or contaminant mitigation techniques. Fractal dimension is useful for predicting the travel distance of fluid particles between two points and the pump time necessary for extraction of fluid, or of a volume of contaminated water, in pump and treat methods (Wheatcraft and Tyler, 1988).

These statistical parameters may be useful for oil and gas field development in the petroleum industry. Variogram interpretation of 3-d seismic surveys could yield contour maps of percent of horizontal variation of a producing formation from a known point to assist in risked reserve estimation. Additionally, pump time to remove a volume of oil or gas could be predicted.

This study differs from previous work in that we calculate the variogram from the seismic section and use it for the interpretation, where power spectra have been the standard tools to date for statistical seismic analysis. For the shorter station intervals utilized in high-resolution seismic surveys for environmental and engineering applications, variograms are more easily interpreted. This study being based on such a survey, we have relied primarily on the variograms for interpretation. However, the use of both methods for extraction of these parameters does enhance the interpretation by adding confirmation of the results.


METHODS

Data collection

We designed and acquired a nontraditional, high-resolution seismic survey to image the near-vertical face of a bench of an open-pit quarry, as a vertically oriented target reflector. We chose this target so future permeameter studies of the exposed face may verify the results of this study. Some 3-d coverage was achieved at a very low cost.

A Bison 9000 12-channel recording system was used for recording. A 1.4 kg single-jack hammer striking a 12-cm hardened steel milling ball was used as a source, with 10 summed strikes per record. Recording twelve foreshot and twelve backshot records at offsets varying from 11.0 m to 1.8 m for each geophone setup resulted in nominally 24-fold coverage. In total, we acquired four lines; one for a length of 43.9 m parallel to the exposed face of the pit at a distance of approximately 15.2 m from the face; and three lines transverse to the face and the first line, tying to the first line and ending at the face (fig. 1).

The line parallel to the face utilized a 0.61-m source spacing. Geophone setups were recorded with foreshots and backshots, then "moved along" six stations and the process repeated. The resulting 24-fold seismic line has a CMP spacing of 0.31 m. Acquisition of the transverse lines utilized a 0.31 m geophone spacing and a 0.61 m source spacing. The result is 24-fold coverage with a CMP spacing of 0.15 m.


Processing

We first generated a suite of 19 constant velocity (CV) stacks from each line, for velocity analysis, giving special attention to low velocities (350 to 550 m/s at intervals of 50 m/s). We expected the velocity of the high porosity, unsaturated diatomite to be in this range. As a result of this analysis, we used a velocity of 400 m/s for migration in the zone of interest. CV stacks generated with velocities in a medium velocity range (nine stacks in the 700-2300 m/s range) helped identify the base of the diatomite at 30 m depth. Additional high-velocity brute stacks (five stacks in the 3000-11000 m/s range) helped identify the "sideswipe" arrivals originating from the vertical face. As was expected these higher-velocity stacks enhanced the sideswipe reflector of interest. Additionally, a single-fold section with increasing offset was generated along the trace of line B (fig. 2). This appears similar to a VSP between Line A and the vertical-face reflector and provides excellent velocity control and identification of the sideswipe event from the face.

We applied both two-dimensional and three-dimensional prestack migrations to the data set. For three-dimensional migration we used the Kirchhoff-summation method of Louie et al. (1988). This method avoids the strong "dip filtering" effect of standard CMP stacking on steeply dipping events such as the face of interest for this study. However, we did not apply any obliquity factor, operator anti-aliasing, or data weighting based on acquisition geometry (Gray et al., 1999). Our migrations thus allow imaging of near-vertical and off-line structures, such as the quarry face, at the expense of accurately reconstructing vertical versus horizontal variations.

Using this migration, we imaged a section in a position along the location of the face (fig.3) and imaged sections along the lines of acquisition of data (such as figure 4). After imaging the face to its proper position, features can be seen and compared between the 3-d migrated face section (fig. 3) and a digital photograph of the face (fig. 5). This migrated section together with those imaged in the transverse direction were used as input for calculating the variograms and power spectra.

Recent work with deep crustal reflectors has resulted in extraction of correlation length and fractal dimension (Pullammanappallil et al., 1996). A relationship between the horizontal reflectivity correlation length and the correlation length of subsurface velocity variations has also been established in recent studies (Levander et al., 1994). A relationship between seismic reflectivity and hydraulic conductivity on a field scale has been established from the literature (Mela, 1997) following work by Gelhar et al.(1979), Lindseth (1979), Schlumberger (1972a,b), and Wylie and Rose (1950). This work shows that the stochastic parameters would be similar for the properties of seismic reflectivity and hydraulic conductivity.


Extraction of variograms and power spectra

A variogram is a geostatistical method of comparing similarity of a data value to neighboring values within a field of data. The variogram is calculated by:


1

(where N is the number of neighboring data points within the lag distance specified, z(xi) is the physical property parameter value of the initial point, and z(xi+h) is the parameter value of the neighboring point being compared).


This value is then plotted against the lag distance between the initial point and the compared points. The resulting plot results in a curve where the variogram value increases with distance to a maximum, and levels off at a lag distance where the total variability of the data field is reached. The geostatistical term for this distance is the range and is the correlation length discussed here.

In addition to this "classical" variogram, we plot a moving-window variogram estimator. Inclusion of a number of values from our band-limited data within a window, instead of a datum from a single point, the moving window estimator developed by Li and Lake (1994) results in a smoother, more stable plot that is more easily interpreted. We plotted both the classical variogram and moving window estimator in the vertical (z) and horizontal (x or y) directions for all migrated sections. These linear plots can be used to interpret correlation length. Carr (1995) demonstrates how log-log plots of such data can be used to determine fractal dimension.

We extracted power spectra from all depth-migrated sections using fast Fourier transform methods (see Claerbout, 1992) in both the horizontal and vertical directions for these sections. This method results in peaks at the spatial frequencies corresponding to the highest correlation. The slopes of such power spectra can also determine fractal dimension. Here, for ease of interpretation, the spatial frequency has been converted to the corresponding wavelength using the diatomite's velocity of 400 m/s. In this way, variograms and power spectra can be used to graphically estimate correlation length and fractal dimension.


STATISTICAL ANALYSIS


The quarry bench under study consists of thin, horizontally bedded (several mm to several cm), fractured diatomite with a few very thin (on the order of several mm) silt and fine-grained sand beds. These diatomite beds horizontally overlay a volcanic ash and tuff that is well below the volume of interest. The result is that the most coherent reflection in the dataset is the bench-face free-surface reflector, which is the emphasis of this study. The vertical components of the section power spectra and variograms show a short range (or correlation length) masked by the variation from the seismic signal. The horizontal component shows this short range originating in the seismic source together with a considerably longer range responding to horizontal lithologic variation.

Figure 6 is a plot of both the classical and moving window estimators of the variogram for the face line in the horizontal direction. Figure 7 shows the vertical variograms, presented here only to demonstrate the seismic signal contribution, the vertical variation being of little interest for the horizontal movement of fluids. The single, steep slope peaking at 3 m lag distance seen on the vertical plot (fig. 7) is characteristic of our band-limited seismic source signal, having a principal wavelength of 3-5 m (as seen on figs. 2 and 10).

The horizontal component of the variogram (fig. 6) contributes most to the interpretation. This variogram shows good development of two distinct slopes before reaching a sill. This is indicative of variance contribution from two sources, the seismic source signal and horizontal lithologic variation. The first steep slope and change in slope at 3 m lag distance corresponds well to those in the vertical variogram. The second, gentler slope from a lag distance of 3 m to 18 m is the contribution of lithologic or reflectivity variation to the variogram. This second slope at larger lag distances is useful for interpreting a correlation length and fractal dimension for the Tertiary diatomite lithology sampled. The vertical variogram (fig. 7) is flat at lag distances greater than 3 m, not showing this second slope. It shows variation in the range expected for the seismic signal and the thin bedding of the diatomite and silt layers, with no longer-length variation.

Correlation length can be visually interpreted from the horizontal variogram by picking the lag where the sill is reached. This lag distance (or range) is the correlation length of the sampled stratigraphic sequence. A point corresponding to 100 percent variance (point B on fig. 6) could be interpreted as correlation length. This point perhaps would be at too large a lag distance, given the long asymptote exhibited here. A point corresponding to the change in slope to the more gradual asymptote, as seen on all the horizontal variograms (point A on fig. 6), could also be used for interpretation, but perhaps would yield a shorter lag distance than desired. We use a lag halfway between these points (designated by the asterisk) for interpretation here. The correlation length estimated using this procedure is 18 m.

A method of extracting fractal dimension (Carr, 1995) utilizes the slope of a log-log plot of the variogram from the relationship:

D = 2 - S/2 2)

where S is the slope of the variogram from this plot. Figure 8 is a plot of the variogram data on a log-log scale. After deleting the first 3.1 m of the lag length (the portion of the variogram where the variation is dominated by the seismic source signal) we can more easily interpret the slope from the lithologic or reflectivity contribution of this plot. Fitting a slope to this plot by eye, the horizontal fractal dimension taken from the face line is 1.94. We have included comparison slopes that would relate to selected fractal dimensions in the lower left corner of figures 10 and 11 to help with interpretation. We obtained similar results (not shown here) for the horizontal variograms calculated from line A for both the correlation length and fractal dimension. These results represent the northing values for these parameters.

Variogram shape is similar for the transverse lines (lines B, C, and D) within the limits of the data collection operation. The longest of these lines (line D at 21 m, not shown) displays two distinct slopes, similar to the face line. The correlation length for all the transverse lines is slightly shorter at an average of 13 m and the fractal dimension averages 1.86, but we consider this good agreement among these lines. These values are the easting component of the parameters.

Figure 9 is the horizontal component of the variogram calculated from the portion of the digital picture containing just the face from fig. 5. The statistics of the photograph depend on the pattern of shadows and color variations on the face. Thus, bulk strength variations and gross layering that the seismic data are also sensitive to dominate the statistics. Using the same methods of interpretation as used for figure 6 results in a horizontal correlation length of 17 m and a fractal dimension of 1.77. While the correlation length compares well, the fractal dimension of the picture is noticeably lower than that of the seismic data. This could well be due to background noise contained in the seismic data that is not present in the digital picture.

The more conventional power-spectral analysis helps elucidate the differences between the horizontal and vertical components. The horizontal and vertical power spectra extracted from the face line are plotted together in figure 10. The contribution of the seismic source wave to the vertical component can readily be identified. Note that the shorter wavelengths (< 8 m) of the seismic signal were not completely imaged by the simple migration of surface data, and dominate the vertical component. This 3 m wavelength is clear in the seismic sections of figures 3 and 4. The 3 m wavelength of the seismic signal corresponds to point A on the power spectra (figure 10) which is the high-frequency contribution of the seismic signal. We applied a low cut filter of 64 Hz during recording (relating to a wavelength of approximately 6 m) and used a 2000 Hz high cutoff. The vertical component shows that the major contribution of the seismic source is at less than 175 Hz (2 m wavelength and longer).

Fractal dimension can be interpreted from power spectra by first determining the slope of the spectrum. Carr (1995) uses:

2


where D is the fractal dimension and b is the absolute value of the slope measured from the power spectrum.

At the longer wavelengths, the horizontal and vertical components of the power spectra separate. The horizontal component remains relatively flat (fractal dimension of 2.5) when wavelengths exceed 8 m. This contribution to the data can only be due to lithologic and reflectivity properties. It is this range of wavelengths that we compare against the variogram interpretation. Figure 11 shows the power spectra extracted from the picture of the face (fig.5). Note the isolated peak of figure 10 (in the range of 2 m to 12 m for the vertical component) is not present in figure 11. The peak is the contribution of the incompletely migrated seismic signal to the power spectrum.

Correlation length can be interpreted from power spectra by determining the length where correlation falls off from the peak corresponding to the lithologic and reflectivity properties. This peak can be seen at approximately 20 m. The correlation length can therefore be interpreted as between 20 m and 40 m, as the next point of control on the power spectra greater than 20 m is 40 m.

While we did not use power spectra for interpreting fractal dimension directly here, a line corresponding to the fractal dimension interpreted from the face variogram has been placed on the power spectrum of fig. 10. While the slopes of the power spectra in both figures 10 and 11 are difficult to interpret, they are consistent with the results obtained from the variograms. On figures 10 and 11 we include diagrams showing the slopes corresponding to different fractal dimensions, to aid in interpreting these power spectra.

The power spectra also point out how the band-limited nature of the seismic data and images can mask any lithologic variation out of the analysis of the vertical components. The seismic response of the thin bedding planes overlaps the range of the seismic signal. Both the vertical variograms and the spectra show the 3 m principal wavelength of the seismic signal. While the horizontal spectrum of fig. 10 shows a "red" spatial correlation flat to large wavelengths, the vertical spectrum shows a "blue" spatial correlation that slopes negatively and disappears at large wavelengths. The "blue" anti-correlation below the 3 m wavelength appears clearly in the migrated seismic sections (figures 3 and 4) as the principal seismic wavelength. However it is not clear from the prestack migrated sections why the horizontal components should have more long-wavelength correlation than the vertical. The pre-stack migration does not act as a strong "low-dip filter" as a CDP stacked seismic section would (Gibson and Levander, 1988), and figures 3 and 4 clearly contain reflections in all orientations. But, not including complete geometrical weighting, anti-aliasing, or obliquity (Gray et al., 1999) in our imaging procedure may have impeded our ability to examine vertical statistics.

Our observations cannot explain why the long-wavelength vertical and horizontal seismic spectra do not agree, while the spectra from the photograph do agree. We can only suggest that the photograph directly samples variation across the face in both horizontal and vertical directions. The seismic data, on the other hand, were recorded by receivers placed at many locations horizontally, but only one level vertically. Thus, the seismic data directly sampled long horizontal wavelengths, with sufficient accuracy that we could recover shades of variation from variograms and spectra. Without a VSP, however, we had only an indirect sampling of vertical variation, through the seismic waves, with noise possibly overpowering the subtle long wavelengths of the vertical variations.


DISCUSSION AND CONCLUSIONS

Comparison of methods

The horizontal correlation lengths of stratigraphic zones can be interpreted by use of variograms as presented above. The spherical model variogram (Carr, 1995), in which the variation rises along a single slope until a maximum variation is reached for the volume under study, best fits the seismic images analyzed in this study as well as most geologic data (Carr et al., 1984). While we used a correlation length corresponding to nearly 100 percent variation for this study, any arbitrary percentage could be used. We could therefore determine lag distances corresponding to any proportion of total horizontal variation. We could then contour these lag distances resulting in a contour map corresponding to percent of horizontal variation from any known point (i.e. a producing well) within a seismic data volume.

We did not use power spectra for extracting a correlation length here, but the spectra are useful in differentiating between lithologic variation and seismic signal. Use of these plots helps to confirm the variation due to lithology and seismic signal. Both appear in the power spectra and can be interpreted from our sections.

Although we can extract fractal dimension estimates from moving-window variograms of seismic images (figure 8), that can be interpreted by the spectra if not visible in them (figures 10 and 11), the accuracy of these estimates is not clear. The images of figures 3 and 4 show the distribution of the seismic reflectivity property, but are severely band-limited by the seismic source, transmission, recording and migration processes. Since the reflectivity property is related to the spatial gradient of rock properties such as porosity and seismic velocity, the spectrum of perfect broadband seismic data ought to be a "blue" or high-pass filtered (f1) version of the rock-property spectrum. Some of this "bluing" is canceled by the low-pass filtering inherent in our simple Kirchhoff-sum migration, since it is an integral transform (Le Bras and Clayton, 1987). Our correlation length estimates are not subject to this "bluing" problem.

In fact, our well-sampled horizontal spectrum at low frequencies is flat (figure 10), whereas the spectrum of the photo has a positive slope representing a high-pass shift from "red" to "blue". Our fractal dimension estimates from seismic reflectivity are also higher, at 1.86 - 1.94, than the 1.77 estimated from the photograph, which also shows spectral "bluing". Although we show here that we can use a moving-window variogram to easily estimate the fractal dimension of reflectivity, additional experiments will be needed before the correlation between reflectivity and porosity or velocity can be established.

With these qualifications, both power spectra and variograms can be used to extract fractal dimension. This statistic can then be used for estimating the pump time for removing a volume of liquid in pump and treat methods for contaminated aquifers (Wheatcraft and Tyler, 1988). The useful distance for analysis of our dataset falls in the range of 6 m to 45 m (low spatial cutoff frequency of the recording instruments and the full length of the line). With values calculated from spatial frequencies at multiples of 0.0128/meter, the resulting power spectra displayed in the range of interest for this study are controlled by only nine values for the face line. Doubling the length of the line to 90 m would result in only one added spectral point, certainly a drawback for the spectral method of determining fractal dimension. Where only nine values contribute to the power spectra analysis used for fractal dimension, 131 values contribute to the analysis of fractal dimension with the variogram method.

A well-developed slope can be seen on the variograms of four of the five seismic lines utilized for this study. The only line difficult to determine the change in slope is one in which the correlation length and the length of the line are nearly the same. Observing the strong agreement of the fractal dimension interpreted from different lines in this manner boosts confidence in this method. While fractal dimensions of similar value are consistent with the power spectra, the scattered values of the sparse data in the zone of the frequency plot utilized makes interpretation difficult. This considerable variation in the slopes of the power spectra results in fractal dimension interpretations from power spectra of questionable value, without first interpreting fractal dimension from a variogram.



Summary

Power spectra and variograms extracted from seismic sections can be useful tools for interpreting the geostatistical parameters of horizontal correlation length and fractal dimension. Extracted from seismic attributes, these parameters can be related to lithologic variations. Geostatistical tools can be applied in such areas of study as hydrogeologic modeling, or reservoir analysis in the petroleum industry. While power spectra and variograms differ in their method of extracting frequency and wavelength information, these differences can serve to enhance interpretation.

Previous applications of extracting seismic attributes have included interpretation of rock type, porosity, and even fluid type. While interpretation of permeability on a broad scale from seismic data has not proven reliable, on a local field scale within a single rock type interpretation of permeability (and therefore hydraulic conductivity) and its stochastic parameters shows promise (Mela, 1997). Using these statistical approaches could lead to future contour mapping of horizontal permeability variation (or hydraulic conductivity) from a known point such as a well (either a producing oil well or a water well). Continued research in this area will benefit hydrogeologic and petroleum reservoir studies.

REFERENCES




ACKNOWLEDGEMENTS

The authors wish to acknowledge the Eagle Picher Corp. and their operations manager, Myron Burdette, who provided access to a site with a geometry that fit the needs of this study so well. Ken Taylor of the Desert Research Institute provided use of seismic instruments for recording the data. Rasool Anooshehpoor provided a computer for field use. Larry Lake made his moving window semivariogram estimator program available, which allowed for easier interpretation f the parameters under study. Dave Benson helped modify this program for this specific use. Finally our data acquisition crew, Nancy Alvarez, Pat Leeper, Matt Mela, and Wuyang Um provided the dedicated work necessary to acquire the high quality data required of this study.


FIGURE CAPTIONS

Figure 1 Layout of high resolution seismic survey on the mine bench.


Figure 2 Line B raw data record showing negative moveout of the face reflection.


Figure 3 Migrated face line. Note the imaged fracture zone from the middle of the bottom rising to the right, and compare with figure 5.


Figure 4 Migrated transverse line showing reflector from face.


Figure 5 Picture of face used to extract power spectra and variograms. The distance between posts is 40 feet (12.2m). Data was extracted between white lines for analysis.


Figure 6 Horizontal variograms for face line.


Figure 7 Vertical variograms for face line.


Figure 8 Horizontal variogram for face line plotted in log-log.


Figure 9 Horizontal variograms extracted from picture of face.


Figure 10 Power spectra extracted from face line.


Figure 11 Power spectra extracted from picture of face.