Quantifying spatial heterogeneity from images

Visualization techniques are extremely useful for characterizing natural materials with complex spatial structure. Although many powerful imaging modalities exist, simple display of the images often does not convey the underlying spatial structure. Instead, quantitative image analysis can extract the most important features of the imaged object in a manner that is easier to comprehend and to compare from sample to sample. This paper describes the formulation of the heterogeneity spectrum to show the extent of spatial heterogeneity as a function of length scale for all length scales to which a particular measurement is sensitive. This technique is especially relevant for describing materials that simultaneously present spatial heterogeneity at multiple length scales. In this paper, the heterogeneity spectrum is applied for the first time to images from optical microscopy. The spectrum is measured for thin section images of complex carbonate rock cores showing heterogeneity at several length scales in the range 10–10 000 μm.


Introduction
The need to describe spatially heterogeneous materials is ubiquitous in physics, chemistry, geology, materials science and related fields. An important practical example comes from the oil industry: carbonates serve as reservoir rocks for more than half of the world's oil, and the processes of sedimentation and diagenesis cause carbonate rocks to have a complex pore space containing a wide range of pore sizes, microporous grains, and an unusual spatial distribution and connectivity of the pores [1,2]. The presence of spatial heterogeneity at many length scales is a fundamental property of these rocks, and this heterogeneity dictates fluid transport in carbonate oil reservoirs. It is natural to use imaging techniques to characterize such spatial heterogeneity: microscopy, x-ray tomography and magnetic resonance imaging are common [3]- [8]. Such images have been used for numerical simulation of fluid flow as well as for constructing statistical models of the pore space [9]- [11].
This paper discusses a different approach to measure the statistical structure in spatial images. The motivation is that although an image carries a tremendous amount of information, it is difficult for an observer to comprehend and quantify this information and to compare different images. Two-point geostatistics represents a classic means of describing spatial heterogeneity quantitatively through the use of tools like variograms and correlograms; however, common geostatistical techniques are ill-suited to describe materials in which a wide range of length scales of heterogeneity occur simultaneously [12,13]. Here, we present a method that expands upon conventional geostatistics to analyze the structure of carbonate rocks displaying heterogeneity at multiple length scales. This treatment rigorously takes into consideration the limited experimental spatial resolution and sample size. Both effects are important in interpreting the images because, intuitively, the measurement's sensitivity is low for heterogeneity at length scales smaller than the resolution and larger than the sample size.
The output of this analysis is the heterogeneity spectrum, which shows the extent of spatial heterogeneity as a function of length scale covering all length scales to which the measurement is sensitive [14,15]. The range of length scales to which the measurement is sensitive is determined quantitatively from the resolution and sample size. The heterogeneity spectrum provides a more comprehensive characterization of naturally occurring media than is accessible from traditional interpretations of two-point geostatistics. Resolving the extent of heterogeneity at different length scales is especially informative for materials that simultaneously present heterogeneity at multiple length scales, such as carbonate rocks.
An analogy can be drawn from other spectroscopic techniques. For example, nuclear magnetic resonance (NMR) spectroscopy measures a time-varying free-induction decay signal arising from the procession of many nuclear spins at different frequencies. Interpretations of those data are routinely performed not in the time-domain but rather in the frequency-domain spectrum that results from a Fourier transform of the time-domain data. Even though the timedomain signal and the Fourier-transformed frequency spectrum contain identical information, the frequency spectrum is much easier to understand, because the peaks in the spectrum are correlated with different chemical groups. Similarly, a peak in the heterogeneity spectrum identifies a structural element in the spatial image.
In this paper, we will describe the geostatistical theory showing how a single measurement can reveal the extent of heterogeneity as a function of length scale at multiple length scales. We also show how the sensitivity of a measurement to different length scales of heterogeneity can be determined from two geometrical parameters: the measurement resolution and the sample size accommodated by the measurement device. We will use this method to analyze two-dimensional (2D) images of carbonate rocks from optical microscopy of thin sections.

Thin section imaging
Thin section samples of two carbonate rocks are prepared and imaged by standard techniques. Briefly, cores of 1.5 inch diameter are vacuum impregnated with a blue-dyed epoxy, and one end is polished down to a flat surface. That end is glued to a glass slide, a thin piece of the rock is cut away from the core, and the resulting piece is polished down until that surface is flat and the rock thickness is approximately 30 µm. The thin sections are imaged with a petrographic microscope equipped with a color camera to collect the image. The resolution of the camera combined with the magnification of the microscope leads to a pixel size of 20 µm. The images are then binarized, with all pixels surpassing a threshold amount of blue assigned as pore (porosity = 100 porosity units (pu)) and all other pixels assigned as grain (porosity = 0 pu). Figure 1 shows the resulting binarized images of the two samples.
Thin section images allow facile observation of relatively large pores. However, it is well known that pore throats and pore bodies much smaller than 30 µm cannot be detected with this technique; higher resolution imaging techniques, such as confocal microscopy, are useful when information on these smaller pores is desired. The goal of this study is not to overcome these general limitations of thin sections. Rather, we attempt to use the images to assess the relevant length scales of spatial heterogeneity in the samples. Therefore, the thin section images are treated as a spatial map showing the location of the larger pores, and the calculated heterogeneity spectrum determines the length scales of heterogeneity relevant to these larger pores.

Inversion method
The process of inverting a spatially resolved measurement into a heterogeneity spectrum has been described previously [14,15]. Here, we will summarize the method and we will detail only the aspects that are unique to this application. The initial step in computing the heterogeneity spectrum involves computing an experimental variogram from a spatially resolved measurement. For a spatially resolved measurement of a given quantity, the experimental variogram is defined as: where N represents the number of pairs of measurements, h is the lag separating the points in the pair, z(x) is the quantity z measured at position x, and the sum runs over all pairs with points separated by h. Conventional methods often simply fit this experimental variogram to a model function such as the exponential or spherical functions in order to obtain the single dominant length scale [12,13]. However, since our samples are known to exhibit multiple length scales of heterogeneity, such a simple interpretation is insufficient. This multiple length scale heterogeneity can be observed qualitatively by simple visual inspection of the images (figure 1), in which a wide range of pore sizes is apparent. The lag is a vector quantity, and variograms can, in principle, be computed along different directions. For the two samples studies here, variograms computed along two orthogonal directions were found to be nearly identical. Such a result typically indicates an isotropic sample, and the isotropic variogram is computed by averaging the two anisotropic variograms; that procedure was followed here, and the resulting variograms are shown in figure 2 [12,13]. In this figure, the variograms are presented on a semi-log scale, which is necessary because the relatively high resolution and large sample size result in an experimental variogram covering almost three orders of magnitude in length scale. The variogram is computed out to 10000 µm, approximately half the sample size, in order to omit points on the variogram containing a small number of pairs [12].
Computing the heterogeneity spectrum requires expressing the variogram as a linear combination of model functions, with each model function in the sum representing heterogeneity at a single length scale. In the hypothetical case of an experiment with zero noise and infinite resolution (point support), the variogram can be expressed as the following sum: where γ ps (h) is the variogram defined on point support, the individual γ model (h; a i ) are the single length scale point variograms, the a i are the length scales of heterogeneity, and the C i s represent the extent of heterogeneity at each length scale a i . The coefficients C i correspond to the weight of the heterogeneity at the length scale a i , so a plot of C i as a function of a i gives the heterogeneity spectrum of this idealized experiment. For a real measurement with finite resolution, spatial averaging within a voxel (regularization) influences the value of the experimental variogram, principally by averaging out heterogeneity at length scales smaller than the voxel size. The model variograms can be regularized analytically according to: In this equation,γ model corresponds to the model variogram averaged over the voxel volume, and regions v and v h are separated by h [12]. The inversion process is general to this point, but for a specific application, a particular model function must be chosen. For application to the thin section images, the common exponential model is selected: From this model function, the regularized model functions are computed according to equations (3) and (4). In 1D, these integrals can be solved analytically. For h l, where l is the length of the 1D pixel, the result is Otherwise, γ model (h = 0; a) = 0 and γ model, reg (h; a) cannot be measured in the range 0 < h < l, which would correspond to separations smaller than the pixel size. Just as the variogram defined on point support was written as a linear combination of model functions, the real experimental variogram can be written as a linear combination of the regularized model functions. For a variogram with discrete points at a finite number of lags and model functions evaluated at a finite value of length scales, the experimental variogram can be related to the heterogeneity spectrum and the γ model, reg according to the following matrix equation: 6 where γ expt (h) is the regularized experimental variogram and is an n-element vector of measured points, γ model, reg (h; a) is an m-by-n matrix relating the length scales of heterogeneity to points on the measured variogram as determined from the integrals listed above, and C(a) is an m-element vector representing the heterogeneity spectrum. After measuring γ expt (h) and computing γ model, reg (h; a), equation (7) can be solved for C(a) using standard matrix inversion techniques, yielding the heterogeneity spectrum. In this inversion, the γ model, reg (h; a) act as basis functions. These functions are the key component in inverting the variograms to the heterogeneity spectra and are shown in figure 3. As basis functions, these functions represent the sensitivity of the experimental variogram to heterogeneity at different length scales. Consider γ model, reg (h = 20 µm; a), which represents the sensitivity of the h = 20 µm point on the variogram (mean squared difference of all pairs of adjacent voxels). If the rock contains heterogeneity at the 1 µm length scale, that heterogeneity will be almost completely averaged out over a pixel, so even the variogram point at this small lag is not sensitive to this length scale of heterogeneity. At the other extreme, if the rock contains heterogeneity at the 100 µm length scale, that heterogeneity will lead to hardly any variation among adjacent voxels, so the variogram point at this small lag is not sensitive to this length scale of heterogeneity, either. In contrast, heterogeneity at the 10 µm length scale is not strongly averaged out over a pixel and does lead to significant variance for adjacent pixels, so this point on the variogram is sensitive to heterogeneity at that length scale. Considering instead γ model, reg (h = 10 000 µm; a), it is found that this point is sensitive to heterogeneity at all length scales that can be observed at shorter lag. Those length scales of heterogeneity will have reached their asymptote at this large lag, so they can be observed with limitation only due to averaging within a pixel. The longest length scale observable is determined as above. Hence, the range of length scales to which the experiment is sensitive is determined on the short end by the size of the pixel (i.e. the resolution) and the long end by the longest lag at which a point on the variogram can be determined (i.e. the sample size).
It is clear from figure 3 that there is considerable overlap in the γ model, reg (h; a) functions. The fact that these functions are not orthonormal presents difficulties in inverting the experimental variogram to the heterogeneity spectrum. Other model functions could be selected, but this set of exponential models was selected for physical reasons. In practice, variograms for these and similar materials typically change slowly. Conceptually, this slowly changing variogram is considered to result from a single characteristic length scale of heterogeneity. The usefulness of exponential and spherical models for fitting experimental variograms results from this simple conceptual understanding. Choosing exponential models here retains this physical interpretation of length scales of heterogeneity, allowing the interpretation of the heterogeneity spectrum as the extent of spatial heterogeneity at each of these characteristic length scales; this choice sacrifices mathematical convenience in the inversion because these functions are not orthogonal. This is similar to NMR relaxometry, where the time constant for the exponential decay signal is indicative of molecular species or local environment. Additionally, the functions are not normalized, reflecting the varying sensitivity to heterogeneity at different length scales.

Structural analysis
Experimental variograms for these two samples are inverted to yield the heterogeneity spectra shown in figure 4. These spectra quantitatively present the extent of spatial heterogeneity as a function of length scale for the larger pore bodies in these carbonate samples. The thin section images contain information over a large range of length scales because adjacent voxels are separated by only 20 µm (center-to-center), whereas voxels on opposite sides of the image are separated by over 30000 µm. Hence, the image is sensitive to heterogeneity at length scales ranging over four orders of magnitude. More quantitatively, the range of length scales of heterogeneity to which the measurement is sensitive can be determined from the sensitivity curves shown in figure 3. Points on the experimental variogram with lag near 10000 µm are sensitive to at least 15% of the spatial heterogeneity at all length scales from 2 to 63000 µm. Assuming that the experimental signal-to-noise ratio allows for the detection of heterogeneity at length scales where the sensitivity exceeds 15%, the range 2-63000 µm was selected as the range of length scales observable in the thin section images. The range includes length scales both below the voxel size and above the sample size because a smoothly varying exponential function was selected to model the heterogeneity; if a square wave model had been selected instead of an exponential model, points on the heterogeneity spectrum below the pixel size and above the sample size would not have been observable (but the physical interpretation of the heterogeneity spectrum in terms of length scales of heterogeneity would have been lost).
Inverting the variograms into heterogeneity spectra also requires selecting the number of length scales (discrete values of a) that can be inverted for simultaneously. Because of the overlap between the sensitivity curves at different lags (figure 3), using an excessive number of length scales can make the inversion numerically unstable. This is especially problematic in the case of an oscillating variogram, which may result from measurement noise or from periodicity in the sample. Reducing the number of length scales helps make the inversion more robust at the expense of resolution in the heterogeneity spectrum. Hence, a compromise must be made in which a sufficient number of length scales is selected such that the heterogeneity spectrum has good resolution without using too many length scales, which would result in an unstable inversion and a heterogeneity spectrum showing great differences in extent of heterogeneity for similar length scales. For these images, two length scales per order of magnitude was selected as a good compromise.
The heterogeneity spectra presented in figure 4 result from inverting the variogram using two length scales of heterogeneity per order of magnitude over the range 2-63000 µm. As previously mentioned, a limitation of thin section images is the difficulty in observing relatively small pores. Indeed, heterogeneity at short length scales in general is difficult to detect in these images because of the thresholding process. Assigning each pixel as either pure pore or pure grain effectively discards all heterogeneity within a pixel. This effect is demonstrated in the heterogeneity spectra, which show hardly any heterogeneity below the pixel dimension. The small values in the heterogeneity spectra at these length scales are not necessarily due to the lack of heterogeneity over this range but rather a result of discarding that heterogeneity during the image processing. If an imaging modality that did not involve thresholding were employed, heterogeneity at a length scale lower than the pixel size by a small amount could be observed in the heterogeneity spectrum; magnetic resonance imaging provides a good example of such an imaging modality [14,15].
The heterogeneity spectra show quantitatively the amount of heterogeneity as a function of length scale. The spectra are shown in absolute units, and perhaps the most striking difference between these two samples is the greater overall heterogeneity found in sample 77. Traditionally, the overall heterogeneity corresponds to the sill (long lag asymptote) of the variogram; however, the sill of the variogram also depends on the experimental resolution because the variogram is not equally sensitive to heterogeneity at all length scales (see figure 3). A better measure of the overall heterogeneity comes from integrating the heterogeneity spectrum, as the variation in sensitivity is corrected. By resolving the contribution of different length scales of heterogeneity, the heterogeneity spectrum contains much more information than the overall heterogeneity. For both samples presented here, the heterogeneity spectra indicate that the most prominent length scale of heterogeneity is approximately 60 µm. For both samples, heterogeneity at length scales in the range 10-1000 µm is relevant, although the spectrum of sample 83 is smoother than the spectrum of sample 77.

Conclusions
The need to visualize spatial characteristics is ubiquitous in studies of materials. The experimental detection of such spatial heterogeneity is always limited by the resolution and sample volume of the specific technique. We develop two-point geostatistics to create the heterogeneity spectrum, which quantifies the extent of spatial heterogeneity as a function of length scale. Key to this method is developing an understanding of the variation of experimental sensitivity to heterogeneity at different length scales. We apply this method to optical thin-section images of rock samples to demonstrate the quantification of heterogeneity for this important visualization technique. This method may also be applied to other imaging techniques and/or images of other samples presenting heterogeneity at multiple length scales simultaneously.