Applying Geostatistics for Investigations of Forest Ecosystems Using Remote Sensing Imagery

Geostatistically based methods that utilize textural information are frequently used to analyze remote sensing (RS) images. The role of these methods in analyzing forested areas increased rapidly during the last several years following advancements in highresolution RS technology. The results of numerous applications of geostatistical methods for processing RS forest images are encouraging. This paper summarizes such results. Three closely related topics are reviewed: 1) specific properties of geostatistical measures of spatial variability calculated from digital images of forested areas, 2) determination of biophysical forest parameters using semivariograms and characterization of forest ecosystem structure at the stand level, and 3) forest classification methods based on spatial information.


Introduction
Forests cover approximately 40% of the global land surface (Westoby 1989) and are one of the most important ecosystems on earth.Because of their economic and environmental importance forests are the subject of intensive monitoring and studies.Field measurements are one of the most exact ways to collect detailed information about forest attributes.However, these measurements are expensive, time consuming and impractical for large areas.Satellite remote sensing (RS) is an efficient and relatively inexpensive way to perform forest studies in a timely and consistent manner over large areas.This method allows acquisition of accurate data that can be easily integrated with geographical information systems (GIS).RS technique is the most suitable for large-scale forest inventories.It may also give important forest and ecosystem information, especially to assess biophysical parameters such as leaf area index (LAI), biomass or net primary productivity.
An important limitation of RS measurements is that data are usually collected at a single spatial resolution, appropriate, for instance, for general land cover observation and analysis, while the forest structure possesses many different scales of variation connected with specific objects.To obtain more detailed results using RS on forest ecosystems, higher spatial resolutions are necessary.Aerial photography can provide data at a variety of spatial resolutions, but such measurements are expensive.The expected availability of high-resolution or multi-resolution RS data in the near future (e.g.Abiodun 1998;Aplin et al. 1999;Atkinson 1999) undoubtedly will increase the role of RS in forest investigations.
Studies involving RS data are based on analysis of both spectral and textural properties of RS images.The use of spectral information has been studied extensively and is the primary basis for estimating properties of the forest and its structure.Up to now, because of the relatively low resolution of RS instruments, texture has been utilized in RS as ancillary information.However, information at pixels of RS images is usually not independent of their neighbors, but rather spatially correlated with them.This correlation may be quantified and used within a classification scheme.
The progress in high-resolution RS makes the role of textural information more and more important.The growing number of high ground resolution satellite sensors has made spatial information on forest properties less expensive and more accurate.The spatial resolution of small satellites has become comparable with single tree canopy diameters.This allows for exact recognition of trees and broad vegetative types as well as detailed characterization of forest ecosystems and stand types at a local level.Textural information is essential when spectral characteristics of different forest stands are similar.The same is true when spectral characteristics of individual stands are highly variable.In such situations, it is advisable to combine spectral and textural information from an image to characterize or classify forest stands.
Many techniques for the investigation of RS images have already been developed to take advantage of texture-based information (e.g., Haralick and Shapiro 1992;Schaale et al. 2001).Of particular importance are geostatistical methods which provide the most flexible approach to spatial correlation analysis.This paper summarizes the use of geostatistical methods for RS studies of forest ecosystems.Because the texture of the forests can be highly complex and can exhibit higher variability than other RS entities, such as geological deposits, sea temperature, etc., the results from the application of geostatistical models for RS of forest communities are often more complicated than for other phenomena.Three closely related topics are reviewed in this paper: 1) specific properties of geostatistical measures of spatial variability calculated from RS images of forest areas, 2) determination of biophysical forest parameters using semivariograms, and 3) forest classification methods based on spatial information.

Specific Properties of Geostatistical Measures of Spatial Variability
Calculated From RS Images of Forested Areas

Geostatistical Measures of Spatial Variability
Geostatistics comprises many methods for evaluating the autocorrelation of spatial data.The central tool of geostatistics is the semivariance function (semivariogram), which is a measure of spatial continuity.The application of the semivariogram requires that the data meet the intrinsic hypothesis for a regionalized variable (Journel and Huijbregts 1978).This hypothesis states that the expected value of the difference in data is zero for all vectors h separating any two points in the region of interest, and that the covariance in the data is a function only of the vector h between samples (Isaaks and Srivastava 1989).
The experimental semivariance for a vector of separation h is derived by calculating one-half the average squared difference in data values for every pair of data locations separated by h.Semivariances can be calculated separately for different directions or combined to produce an omnidirectional semivariance.These values are then plotted against the distances between data pairs.Such a plot is also commonly referred to as a semivariogram or variogram.The following formula is the most frequently used for the semivariance calculations: where x i is a data location, h is a lag vector, Z(x i ) is the data value at location x i , and N is the number of data pairs spaced a distance and direction h units apart.A common form of the semivariogram is shown in Fig. 1.
Semivariance functions are usually characterized by three parameters: -Sill -the plateau that the semivariogram reaches.
The sill is the sum of total variation explained by the spatial structure and nugget effect (described below).-Range (range of influence or correlation) -the distance at which the semivariogram reaches the sill, or at which two data points are uncorrelated.-Nugget effect -the vertical discontinuity at the origin.The nugget effect is a combination of sampling error and short-scale variation that occurs at a scale smaller than the closest sample spacing.Semivariance calculations can also be performed with data from RS images.The semivariogram is calculated from the raster images using digital numbers (DN) as data values Z(x i ).
A semivariogram with classic form can be an efficient tool in the analysis of RS images.For example, it can be used for estimating appropriate spatial resolution because its range determines the distance above which ground resolution elements are not related (Curran 1988).
Another important measure of spatial correlation is the cross-semivariogram: where x i is a data location, h is a lag vector, Z(x i ) and W(x i ) are the DN values at location x i for different bands, and N is the number of data pairs separated by length of the vector h.The cross-semivariogram quantifies the joint spatial variability (similar to cross-correlation) between two radiometric bands.A similar, but simpler measure of joint spatial variability is the pseudocross-semivariogram: Other measures of spatial variability are described in detail by Deutch and Journel (1998).

Typical Forms of Experimental Semivariograms and Their Interpretation
The classic semivariogram form (shown in Fig. 1) often results when analyzing RS images of forest areas.Such semivariograms are useful for determining a vegetation community structure (Curran 1988).Their range indicates the distance at which pixel properties are no longer correlated, and therefore provides a measure of the size of the elements embedded within the image as expressed by the semivariogram.The range is also often related to the size, or scale, of the largest elements in the scene that produce the correlation structure (Jupp et al. 1988(Jupp et al. , 1989)).It provides a measure of the distance around a point at which spatial interpolation or processing is valid.In the case of RS data, it indicates the optimal spatial resolution for characterizing the elements embedded within the image or optimal window size at which to apply textural measures to the image data (Treitz and Howarth 2000).
The sill provides a measure of the variability of the reflectance values across the RS image of the stand at distances where there is no spatial dependence between the reflectance values.The sill has been associated with the complexity of the image data and, hence, with the complexity of the target surface (i.e., forest canopy).The nugget effect represents spatially independent variability that generally arises from measurement error (Huijbregts 1975), and short-scale variations that occur at a scale smaller than the closest sample spacing (pixel size).
The investigation of RS images of forests requires detailed knowledge of internal and external factors affecting their optical properties.The main external factors are: sun elevation and azimuth, atmospheric scattering, properties of the RS instrument.The main internal factors are: stand type, species of trees, canopy geometry, tree row orientation, distances between trees, underbrush, leaf litter and optical properties of the background.The terrain geometry and the season of the year in which measurements are made are essential for the correct interpretation of results.Guyot et al. (1989) and Civco (1989) discuss this problem in detail.
The semivariograms calculated from RS images of complex forest objects can be much more complicated than the classic ones.They strongly depend on the above-mentioned factors affecting optical properties of images.Semivariances also depend highly on ground spatial resolution of the RS (Atkinson 1993), leading to issues of scale.Decrease of the spatial resolution has the following effect on semivariogram: 1) the height of the sill decreases, 2) the range of influence increases, and 3) the height of the semivariogram at the first lag increases relative to the sill (Woodcock et al. 1988a).Marceau et al. (1994a, b) describe in detail problems of scale in forestry.The scale problem refers to the fact that results computed from the same data depend on the size of the units of analysis, and therefore a variety of different results may be obtained from the same data.
Different forms of semivariogram were observed in practice, such as common (classic), periodic or non-spatial ones (Curran 1988).This first type of semivariogram appears when repetitive patterns are studied, and the second type appears when random patterns are investigated.In the case of spatial resolution from 20 to 30 m, the periodic behavior of the semivariogram is connected mostly with either forest stands or a group of tree properties (Strahler et al. 1986).In the case of spatial resolutions below 5 m, the periodic behavior of the semivariance is connected rather with individual tree observations.There are also often observed "unbounded" forms of semivariograms (especially, when spatial resolution above 20 m was used).The unbounded semivariogram represents a situation in which either a trend or many different-range spatial correlations are observed or correlation is maintained beyond the sampled area.Atypical semivariogram forms have been reported and described in detail by Curran (1988).They are much more difficult to model and interpret.Semivariograms often show many scales of correlations existing in studied forests.In relatively simple cases, such semivariograms can be fitted using a few different models ("nested structures", see Journel and Huijbregts 1978).The ability to derive information about multiple scales of variation in images from semivariograms is one of the most attractive features of semivariograms.On the other hand, the usefulness of complex semivariograms forms for analyses of forest attributes or classification is usually considerably less because of severe difficulties in their parameter estimation and interpretation.

The Effect of Spatial Resolution and Optical Wavelength on Semivariograms with Forest Data
The spatial resolution of RS images is an essential factor that must be taken into account when analyzing semivariograms produced from RS forest data.Spatial resolution refers to the smallest object on the ground that can clearly be "seen" by the sensor (Lillesand and Kiefer 1979).In terms of RS data, one can think of it as the ground area covered by one pixel.According to Strahler et al. (1986), RS data is considered "low-resolution" when the ground area covered by one pixel is larger than the actual objects on the ground.Otherwise, it is considered to be "high-resolution".Problems of spatial resolution and level of plant recognition typical in forest applications were discussed in detail, e.g. by Wulder (1998), Lefsky and Cohen (2003), Tomppo et al. (2003).According to the latter author, in terms of global forest assessments, image resolution is considered to be low if the pixel dimensions are greater than 1000 meters, medium if between 100 m and 1000 m, high if between 10 m and 30 m, and very high if below 5 meters.The high resolution Landsat satellites are particularly useful for environmental research because their spatial resolution (e.g., 30 meters for Landsat TM) usually corresponds well with the scale at which the environmental system of interest operates.This is referred to as the "operational" scale (Sabins 1978).The very highresolution satellites, such as QuickBird (Digital-Globe), OrbView, and IKONOS (ORBIMAGE), having spatial resolutions as low as 61 cm, make possible the evaluation of individual trees.Measurements in RS are not point-specific, but are integrated over some area, referred to as the unit of regularization, which is equal to the field of view of the sensor.The resulting semivariogram is referred to as regularized and is different from the one obtained using point measurements.Usually the resolution -cell size of the RS image -is taken as the unit of regularization.In many cases the effect of regularization on point-based semivariograms can be determined analytically, as described in detail in the geostatistical literature (Webster and Olivier 1990;Isaaks and Srivastava 1989;Journel and Huijbregts 1978;Jupp et al. 1988Jupp et al. , 1989)).However, in the case of RS images, the relationship between the regularized and point semivariogram is much more complicated due to the effect of noise and non-uniform sensor response arising e.g., from geometric distorsion of RS imagery (Atkinson 1993).The problems of regularization of the semivariograms obtained from RS images have been thoroughly discussed in the literature (Woodcock et al. 1988a(Woodcock et al. , 1988b;;Woodcock and Harvard 1992;Atkinson 1993Atkinson , 1997, etc.), etc.).Atkinson examined the behavior of the experimental semivariogram with the spatial resolution as a means of inferring the relation of spatial resolution to information (structural variability) and measurements error (nugget effect).He stated, that measurement error must be identified, removed from the semivariogram, and regularized separately from the underlying variation.
An important consequence of the regularization process is that the semivariogram calculated from RS images within a given cover class is only characteristic of that class for a specified sensor system.Variation at a scale finer than the scale of regularization cannot be detected and variations less than two to three times the scale of regularization cannot be reliably characterized (Woodcock et al. 1988a).
Another fact to remember is that spatial variations of reflectance observed on RS images are strongly dependent on wavelength due to the effects of many phenomena on reflectance of forest areas.Green plants strongly absorb visible electromagnetic radiation and strongly scatter near-infrared radiation (Curran 1980).Therefore, different semivariograms can be obtained observing the same area with constant spatial resolution for different wavelengths.Typically, response of green vegetation is generally higher in the infrared portion of the wavelength spectrum when compared to the visible portion.Therefore, the infrared-based semivariogram gives much more information on forest structure than a semivariogram calculated from the visible light reflectance.One approach is to calculate the semivariograms from vegetation indices such as the Normalized Difference Vegetation Index (NDVI) as in Rouse et al. (1973).The NDVI index has been found to be a sensitive indicator of the presence and condition of green vegetation (Tucker et al. 1986).There are many different vegetation indices that are appropriate for investigating different forest characteristics (see e.g., Elvidge and Chen 1995).However, because of nonlinear dependencies between vegetation indices and biophysical parameters, the usefulness of these indices is sometimes limited (Wulder et al. 1996).For the above-mentioned reasons, the parameters derived from cross-semivariograms that describe multi-band spatial correlations can be interesting alternatives to vegetation indices (Chica-Olmo and Abarca-Hernandez 2000; Jakomulska and Clarke 2001).

Determining Biophysical Forest Parameters and Characterization of Forest Ecosystem Structure at a Stand Level
Forest biophysical parameters are numerical characteristics that enable the calculation of complicated forest properties using a single measure.
The high spatial resolution of new satellites provides greater accuracy in determining biophysical forest parameters.We present here some important recent work making use of geostatistical methods for RS based forest investigation at different spatial resolutions.

Empirical Studies
Analyses of RS images texture of forest stands are often connected with simultaneous ground measurements, especially when a supervised classification is considered.An intensive fieldwork must be conducted to determine parameters of observed stands.Some parameters such as illumination conditions are difficult to control and isolate during analysis.Observations and interpretations of real forest stand features are necessary.Hereafter, we present the most meaningful results.Cohen et al. (1990) used semivariances to exploit spatial information inherent in RS imagery for different coniferous forest stands in the Pacific Northwest region in the USA.The semivariances calculated using DN of digitized aerial videography at a nominal 1 m pixel size were used to evaluate the potential utility of results obtained using SPOT HRV and Landsat TM satellites for canopy structure analysis.The calculations of semivariances were repeated after the video images were spatially degraded to 10 and 30 m pixel sizes, for SPOT HRV and Landsat TM, respectively.The particular objective of this study was to evaluate RS utility for distinguishing among the stands having different canopy structures, which included even-and uneven-aged young, evenand uneven-aged mature, and old growth stands.The authors found that at 1 m spatial resolution the sills corresponded to the presence of vertical layering in the canopies and to percent canopy cover, while the ranges of the semivariograms were closely related to the mean tree canopy sizes of the analyzed stands.Therefore, they concluded that at 1 m resolution the sills of the semivariances are good discriminative parameters for different types of forest stands.Semivariances based on 10 and 30 m pixels contained significantly less useful information.However, calculations by the above-mentioned authors were limited only to the red band of the images and did not combine multiple-band textural information, which could have potentially improved their results.Treitz and Howarth (2000) examined high-resolution data collected at the stand level.Different mono-specific or mixed forest stands using the Compact Airborne Spectrographic Imager (altitude of 600 m, spatial resolution 0.73 by 5.36 m) were studied by means of semivariances.It was concluded that semivariograms calculated from high-resolution data are closely related to some forest parameters.Among other conclusions drawn are the following: a) range estimates differ between various forest ecosystem classes, and appear to be related to the crown diameters; b) range estimates vary as a function of wavelength, i.e., contrasting estimates were derived for the visible and near-infrared wavebands; c) a direct relationship exists between percentage of understory cover and semivariance; d) although the semivariance is related to the density of trees in natural ecosystems, it is more closely related to the whole density of all layers of vegetation.The obtained site-specific results were applied to choosing an optimal spatial resolution for a landscape-scale (1:20 000) forest ecosystem classification.Wallace et al. (2000) obtained similar results from ground measurements of vegetation communities in the Mojave Desert of California.Fitted semivariogram models revealed that different vegetation communities have distinctive spatial properties, according to the model parameters.It was concluded that spatial pattern information produced by ground measurements can improve significantly broad-scale vegetation classifications produced by low-resolution RS systems.Franklin and McDermid (1993) analyzed Digital Compact Airborne Spectrographic Imager (CASI) data and SPOT satellite HRV multi-spectral imagery data in order to establish empirical relationships between RS data and a variety of forest stand and individual tree parameters.Another task was to discriminate between such forest characteristics as volume, density and canopy coverage beyond the scope of traditional medium-scale aerial photo-interpretation by exploiting the specific advantages offered by high-resolution RS imagery.The authors replaced each pixel with the mean of a local window for smoothing.The size of the window was determined using the semivariance.Selecting a proper size both for derivation of texture and for filtering was critical.The window size was customized for specific stands.The authors concluded that analyses of forest stands based on custom windows produced at least a five percent improvement over any combination of static window sizes.Wulder et al. (1998) showed that inclusion of texture, which acts as a surrogate for forest structure, into the relationship between LAI and the normalized difference vegetation index (NDVI) can greatly increase the accuracy of modeled LAI estimates.Compact airborne spectrographic imager (CASI) data were used for this study (at an elevation of 700 meters with 1 m spatial resolution and five user-selected spectral bands).Firstorder and second-order moment texture values (Peddle and Franklin 1991) as well as semivariance moment texture (SMT) values calculated from DN of RS images were examined for their relationship with LAI.SMT values were extracted from unidirectional semivariograms calculated within a moving window.The nugget, sill, range, slope of semivariance, and mean semivariance were used as surrogate measures of image texture due to the characterization of spatial variability in DN.The ability to increase the accuracy of LAI estimates was demonstrated over a range of forest species, densities, closures, and successional regimes.Deciduous stands, spanning an LAI range from about 1.5 to 7, exhibited a moderate initial relationship between LAI and NDVI with a correlation coefficient of 0.42.Inclusion of additional texture statistics into the relationship between LAI and NDVI further increased this correlation coefficient to 0.61, representing an increase in the ability to estimate hardwood forest LAI from RS imagery.Mixed forest stands, which were spectrally diverse, had an insignificant initial correlation coefficient of 0.01 between LAI and NDVI, which was improved to a significant correlation coefficient of 0.44 with the inclusion of the semivariance moment texture.Empirical models for the estimation of LAI using STM at a high resolution may then be used to provide new information to increase the accuracy of estimates at a lower resolution (Wulder et al. 1996;Franklin et al. 1997).For example, the LAI estimates derived from CASI imagery at a high resolution may be used to provide additional information to assess the sub-pixel characteristics of lower resolution data such as Landsat TM.Wallerman et al. (2002) demonstrated a new kriging method (the so-called edge model) for prediction of forest stem volume.This method incorporates spatial information derived by an edge-detection algorithm to increase the prediction accuracy and was applied to Landsat TM data.The basis of the edge model is to use knowledge of present edges to select only the neighboring plots on the same side of the closest edges as the unobserved point to be predicted.The results were compared with global (using the data from the complete forest area) and stratified kriging (using the data from pre-selected strata).Prediction based on the new edge model was superior to prediction using the global kriging and almost as accurate as the stratified prediction method.
Lévesque and King (1999) studied the forest structural damage at an Acid Mine Site in Ontario, Canada, using an airborne digital camera.The principal objective of the study was to check whether semivariance statistics from very highresolution (0.25, 0.5 and 1.0 m) imagery could be used to model forest canopy and individual tree crown structural variations related to various degrees of pollution-induced damage.Sampling was conducted over the forest canopy as well as within individual tree crowns.The semivariance analysis identified several significant relationships between image and forest parameters that are neither evident in single band raw spectral data nor in simple first-order image texture.The authors found that the range and sill of omnidirectional semivariograms were well correlated with a visual tree-stress index.The ranges of isotropic semivariances were correlated with height and were sensitive to tree crown size and canopy closure.The semivariances extracted from 0.5 m pixel images were suitable for mapping structural and textural information related to forest damage at the canopy level, while the semivariograms extracted from 0.25 m pixels depicted information at the tree crown level.
Several other authors have also used semivariograms calculated from RS data to extract information about forest damage.For example Bowers et al. (1994) used semivariogram statistics (range, sill, and nugget) derived from SPOT HRV panchromatic and multispectral data to reveal and model Balsam fir damage caused by the woolly adelgid.

Model-based Studies
The major limitation of approaches based on analyses of RS images of forest areas arises from the fact that it is impossible to control all the factors of forest properties influencing RS images.For this reason, it is not an easy task to interpret RS images and their texture parameters.This is why model-based studies of forests are intensively carried out.The results obtained from models are verified empirically and appropriately calibrated.Woodcock et al. (1988a) undertook a detailed model-based study to identify dependencies between spatial characteristics and modeled images of forests.They used a simple scene model of randomly located overlapping discs on a continuous background.They related directly texture information to parameters that characterize the discs (e.g., size, density, area of background).The authors also studied the influence of different spatial resolutions on spatial statistics.They additionally modeled a coniferous forest scene comprised of four kinds of differently lighted objects: illuminated canopy, shadowed canopy, illuminated background, and shadowed background.Important findings included the following: a) the sill is related to the proportion of the area covered by objects; b) the range is related to the size of the objects in the scene; c) the shape of the semivariogram and the range of influence are more closely related to the area of objects than to their shape; d) the shape of a semivariogram is related to the overall variance in the scene; e) reducing the spatial resolution result in a lower sill and larger range of influence, as well as in faster increase of the semivariogram height at the distance of the first lag.
St-Onge and Cavayas (1995) studied the influence of tree size and density on the spatial structure of high-resolution images of forest stands.As a first step they simulated artificial images with the help of a three-dimensional canopy model and a geometrical ¯optical model (Granberg 1987).Thus, they could analyze the influence of these parameters on image texture and determine prediction equations for crown diameter, tree height, and stand density, based on image texture.They predicted a stand structure by using a set of three functions to relate the directional range parameters to the forest structure parameters.The directional semivariograms proved to be an effective tool for studying structure of the modeled forests.
To check the modeled results they used very highresolution (36 cm) airborne images.Good results were obtained for such attributes as tree height, crown diameter, and stand density.The authors also tested the influence of spatial resolution on semivariograms.They decreased the resolu-tion from 0.36 to 2.16 m by averaging the pixel values in 6 by 6 pixel squares.In this case, good results were obtained for tree height and stand density.The authors concluded that it might be possible to establish an absolute spatial signature of forest stands, which cannot be attained in classical spectral RS.They also concluded that good prediction of forest inventories could be achieved by combining a segmentation algorithm with their method based on the directional semivariograms at a lower spatial resolution of 2 m.
Bruniquel-Pinel and Gastellu-Etchegorry (1998) thoroughly studied the sensitivity of textural information from high resolution RS images of a forest pine plantation in Les Landes, France to a number of biophysical parameters: crown diameter, distance between trees and rows, tree positioning, leaf area index, and tree height.They also investigated the influence of acquisition parameters such as spatial resolution, spectral domain, spectral viewing, and illumination configurations.The authors calculated semivariograms from simulated images generated using the discrete anisotropic radiative transfer (DART) model.This model provides realistic simulations of images of 3-D scenes under any viewing and illumination configurations (Gastellu-Etchegorry et al. 1996).Image texture sensitivity to forest parameters was analyzed in the visible (green) and near-infrared domains, with a "reference parcel" of which biophysical parameters were varied, and optical parameters were kept constant.Results indicated a complex dependence of semivariogram characteristics (range, sill, amplitude of oscillations) on biophysical and acquisition parameters.The authors tested theoretical results against high-resolution airborne data (1.67 m resolution) and they found that the sill and the oscillation amplitude of semivariograms increase strongly as LAI increases in visible light.These two textural parameter values increased considerably for LAI values between 1 and 4 for horizontal variograms (those computed in the direction parallel to tree rows), and between 1 and 5 for vertical variograms (those computed in the direction perpendicular to rows).The influence of LAI was also tested in the case of a so-called "natural forest", a forest where trees are randomly distributed in the scene.Similar to the case of pine stands, sills of the natural forest increased with LAI increases, whereas the range did not show any clear tendency.The sills of the natural forest were larger than those of pine plantations with the same LAI, due to larger radiometric variability.Important differences in texture parameters also existed between simulated and real RS images.These differences were attributed to the atmospheric diffuse scattering and associated environmental effects, RS instrument characteristics, and variability of understory characteristics.

Review of Large-Area Classification Methods Based on Spatial Information
As we already mentioned, the spatial resolution of the RS instrument is essential for the information content of RS images.Low-resolution RS is the most appropriate for large-scale forest classification.American Landsat and French SPOT satellites still provide the majority of RS images in use today.Correct classification methods based on spectral information from Landsat TM images can discriminate between classes having sufficiently different spectral characteristics.For example, it is relatively easy to discriminate between water, field, hardwood and coniferous forests, estimate LAI (e.g., Chen and Cihlar 1996) or observe some forest changes (e.g., Woodcock et al. 2001).Today the 30 m spatial resolution of the Landsat satellite is considered to be the upper limit that is suitable for forest observations (Fazakas and Nilsson 1996;Wulder 1998).Because it is not possible to observe and analyze objects smaller than the spatial resolution of the sensor, the amount of spatial pattern information in semivariograms calculated from low-resolution RS forest images is limited.Landsat TM data combine different spectral contributions within a pixel representing e.g., trees, underbrush and leaf litter.All these contributions can influence semivariances and other spatial measures of correlation.Therefore, studies applying geostatistics in low-resolution RS imagery for forest observations are relatively rare.The earlier empirical study by Cohen et al. (1990) suggested limited applicability of geostatistical methods for Landsat TM forest image analysis.An example of forest inventory com-bining the field data and those from Landsat TM images was given by Lappi (2001).High-resolution RS data are still expensive and time consuming for computer processing.However, the results calculated from high-resolution data give precise information on forest properties for smaller areas.Undoubtedly, the high-resolution orbital sensor will play a significant role in large-area classification in the future.However, this will demand different image processing approaches.Hereafter, we summarize the texturebased classification methods, focusing mainly on the recent achievements in this field.
There are many methods that use statistical and spatial information for classification.Such methods usually fall into three broad categories: -Methods based on non-spatial statistical information -Non-geostatistical methods that exploit spatial autocorrelations -Methods based on geostatistical measures of spatial variability.

Methods Based on Non-Spatial Statistical Information
The most popular methods of this type are based on adding some statistical information as new (ancillary) layers of RS images.Additional layers contain statistical information that is extracted from radiometric bands using appropriate filters.The most common are calculations of the mean, median and variance in a local moving window.The size of the moving windows is usually of a few by a few pixels.The ancillary data become an additional feature of each pixel and can be used for classification purposes (see e.g., Arai 1993;Chen et al. 1997;Haralick and Shapiro 1992;Woodcock and Harvard 1992).

Non-Geostatistical Methods That Exploit Spatial Autocorrelations
There are also many methods based on spatial autocorrelations, which do not employ semivariance analysis.For example, Labovitz and Masuoka (1984) studied the influence of spatial autocorrelations in signature extraction during pixel classification.Goodchild (1986) studied spatial correlations using Moran's "I" and Geary's "c" statistics for analyzing raster images.Anselin (1995)

Methods Based on Geostatistical Measures of Spatial Variability
Below we describe the applicability of geostatistical methods to RS forest classification, focusing on the most recent results.

Methods Based on Pixel by Pixel-Derived Measures of Spatial Continuity
Classification, or clustering, is the process of sorting pixels into a finite number of individual classes or categories of data based on their attributes.If a pixel satisfies a certain set of criteria, it is assigned to the class that corresponds to these criteria (ERDAS Field Guide 1990).
The most recent and effective methods that use textural information are based on computing a set of textural measures of spatial variability (TMSV) using a moving window.Because the choice of an appropriate moving window size is an especially important issue when spatial measures of the image are calculated and modeled, we discuss this problem at the beginning of this section.Most applications presented in the literature use a fixed window size.This is the most straightforward and the fastest technique, but it is not always precise.A trade-off exists between moving window sizes.On the one hand, a window that is too large risks encompassing several classes in one window.On the other hand, a window that is too small does not contain enough information for accurate computation (and eventually modeling) of a semivariogram.Often the moving window size is chosen experimentally during the training phase (e.g., Miranda et al. 1992;Miranda and Carr 1994;Miranda et al. 1996;Chica-Olmo and Abarca-Hernandez 2000).However, the most promising strategy is to determine the moving window size and shape based on the actual image variability.Such a window is often referred to as a geographic window (Merchant 1984;Dillworth et al. 1994).Merchant (1984) proposed using geographic windows for Landsat TM image classification.In order to automate determination of the optimal geographic window size, some authors (e.g., Franklin and McDermid 1993;Franklin et al. 1996) employed semivariances.This type of moving window was also applied in forestry for LAI estimation (Wulder et al. 1998), and recently, for assessing forest landscape structure (Ricotta et al. 2003).
Two different approaches were proposed for combining geostatistical methods with the moving window methodology.The first consists of calculating and accurately modeling the semivariogram in a moving window, and then using some of the modeled semivariogram parameters for class identification.Ramstein and Raffy (1989) were the first to use this approach in analyzing a few distinct surface cover classes (forest, fields, and various urban areas) using isotropic semivariograms of Landsat TM images at a 30 m spatial resolution.They concluded that the sill of a semivariogram is not appropriate for classification and proposed a one-parameter texture classification based on the range of the semivariogram.The experimental semivariograms were fit using an exponential model.In order to reduce computing time, the semivariogram parameters were not obtained from a traditional least-squares procedure, but were calculated from the formulas containing the means of DN and the means of squares of DN as well as the values of the semivariances at the first lag.It was found that the optimum moving window size was 11 by 11 pixels.The classification levels for the range of semivariograms calculated using the Landsat TM red channel were: a < 1.8 for urban areas, 1.8 < a < 3.1 for fields, and a > 3.1 for forests.
In the second and more common geostatis-tical approach to classification semivariograms calculated in moving windows are not modeled.Because of the complicated forms of semivariograms calculated from forest areas, the approaches based on accurate pixel-bypixel modeling of semivariograms in a moving window are not very useful.Large amounts of RS data sources make accurate semivariogram modeling computationally unreasonable.Instead, a set of semivariance-derived measures is used (e.g., semivariance values at consecutive lags, slopes of semivariances for consecutive pairs of lags, mean semivariance, some measures of semivariance shape, etc.).This method was used by Miranda et al. (1992), Miranda andCarr (1994), andMiranda et al. (1996) by employing a semivariance textural classifier (STC) algorithm.In the STC algorithm semivariances for the first six consecutive lags have been used as an additional input layer to image classification.For example, the semivariances were calculated from Shuttle Image Radar-B (SIR-B) data (Band L-23.4 cm wavelength) in order to discriminate vegetation in Borneo (Miranda et al. 1992).Four classes of surface cover type were selected: water, tidal forest, coastal lowland forest and swamp.The classification process consisted of two major steps: training and classification.The training was performed using a 22 by 22 pixel window for each class.During classification, semivariances were calculated within a moving window of 7 by 7 pixels.The results of texture analyses based on classification using the STC algorithm agreed well with the results attained by Ford and Casley (1988) through visual interpretation of the same SIR-B dataset.Similar methods were applied (Miranda et al. 1996) for vegetation discrimination in northwestern Brazil using JERS-1 (Japanese Earth Resources Satellite 1, renamed to Fyuo-1 in 1993) synthetic aperture radar data with an 18 m spatial resolution.The following classes were investigated: water, open vegetation, dense vegetation and flooded vegetation.Classification accuracy was assessed in a confusion matrix built for the training sites.The classification accuracy was 0.58, 0.60, 0.67, and 0.68 for water, flooded vegetation, open vegetation, and dense vegetation, respectively.Carr (1996) presents two computer programs (MXTEXT and MXMULT) that use minimum-distance-to-mean or Bayesian maximum likelihood algorithms for spectral classification, allowing classification of image texture based on the local semivariograms calculated in moving windows.
A similar method for texture-based image classification was developed by Chica-Olmo and Abarca-Hernandez (2000) for RS data classification of two different geological deposits.Although no forests were classified in this study, it seems reasonable to consider this method here.In their study the TMSVs were calculated from the principal components (PC) on the radiometric Landsat TM5 bands.Six Landsat TM bands (all except TM6) were divided into two groups (visible and infrared).The most representative principal components of each group, namely PC1 and PC3, were calculated.Then a madogram and semivariogram as well as a cross-semivariogram and pseudo-cross-semivariogram for the two above-mentioned PC's were calculated in moving windows of 7 by 7 pixels for the first lag only.As a result, a geostatistical texture image of six TMSVs was obtained (semivariograms and madograms for each PC and cross-semivariogram and cross-pseudo-semivariogram between PC1 and PC3.Additionally, variances for each PC and four fractal dimensions along two perpendicular directions for each PC were calculated in a moving window.Different combinations of the above variables were studied to improve classification accuracy with a likelihood decision rule using ERDAS software with the radiometric information (ERDAS Field Guide 1990).Two important conclusions were reached: 1) the textural information improved classification accuracy (by 19.7%), and 2) multivariate measures, such as cross-semivariances, pseudo-cross-semivariances and madograms, played an important role in increasing classification accuracy (by 11.6%).
St-Onge and Cavayas (1997) developed a set of algorithms designed for large-scale inventory of forest structure parameters from high-resolution imagery (≤ 1 m).Texture information was first derived by measuring the range of the semivariogram that was calculated from monochromatic images in three different directions using a moving window.The semivariogram ranges were then used to predict crown diameters, stand density, and crown closure through regression equations on a per-pixel basis.A region-growing algorithm was applied to these three regressionestimated images to identify the limits of the forest stands.Calibration of the prediction equations was made using artificial images created by a geometrical-optical process.It was found that forest stand boundaries can be adequately identified on artificial images and that average forest structure estimates within each delineated stand are close to the actual values.Preliminary use of the proposed method yielded good segmentation and per stand structure estimates when applied to real images acquired with the MEIS-II (Multi-detector Electro-optical Imaging Scanner) airborne sensor.
Jakomulska and Clarke ( 2001) performed a texture-based classification of chaparral vegetation using data from the Airborne Data Acquisition and Registration System 5500 with a 2.5 m ground resolution.They chose green, red and infrared bands to analyze six different cover types.Experimental semivariances, cross-semivariances and pseudo-cross-semivariances between three studied bands were calculated for each image pixel using a technique that incorporated a moving window of a changeable size.The following sets of indices were calculated in a moving window for each pixel image: 1) individual parameters (semivariogram values for different lags); 2) global parameters such as the sum of the semivariances/cross-semivariances, mean semivariogram, sum of the absolute differences between semivariance models and experimental semivariogram, and sum of the absolute differences between semivariograms and cross-semivariograms; 3) parameters connected with the range or sill (difference between semivariogram and cross-semivariogram at a range, ratios of sills for band pairs and range values); and 4) measure of semivariogram shape (e.g., ratio of sill and range).To evaluate the discriminative power of the above-mentioned texture indices, a set of training samples from studied classes were used.Due to the large number of derived parameters, the selection was made through principal components and correlation analyses.It was found that the best discriminative parameters were the sills of the semivariogram and cross-semivariogram, the mean semivariance and cross-semivariance, the sum of absolute differences between spherical semivariograms and semivariograms at lags up to the range, and the pseudo-cross-semivariogram at lag 0 and lag equal to range.The overall accuracy using both spectral and textural information versus only spectral information was 81.7% and 67.5%, respectively.This means that the overall classification accuracy increased by nearly 15% thanks to the inclusion of textural information.It was also found that discrimination between classes found in shaded areas is possible if textural parameters are used.
Cokriging techniques are also promising in forest inventory.Integration of different types of RS information (e.g., multisensor data) can be of even further use.For example, Hudak et al. (2002) integrated the lidar and Landsat data for estimating and mapping forest canopy height in western Oregon.Land detection and ranging (lidar) data provide detailed information on canopy structure in the vertical plane, but over a limited spatial extent (Lefsky et al. 1999).Landsat data provide extensive coverage of generalized forest structural classes in the horizontal plane but are relatively insensitive to variation of forest canopy height.Spatially continuous lidar coverage was sampled in eight systematic patterns to determine which lidar sampling strategy would optimize lidar-Landsat integration.The following five non-spatial and spatial methods were tested for predicting canopy height: regression, kriging, cokriging, and kriging and cokriging of regression residuals.The authors concluded the integrated models that kriged or cokriged regression residuals are the most suitable strategy for estimating and mapping canopy height at locations not sampled by lidar, and that a 250 m discrete point sampling strategy was most efficient for an intensively managed forested landscape in western Oregon.

Approaches for Classification Not Based on Moving Windows Techniques
Some geostatistically based RS image classification methods are not performed on a per-pixel basis.For example, Dungan (1998) studied the applicability of geostatistical methods for spatial prediction of vegetation quantities, such as biomass, canopy cover, stems per area unit, etc., derived from radiance of pixels of RS images using a synthetic example constructed on the basis of RS data.Three geostatistical methods: kriging, cokriging, and stochastic conditional simulation were compared with a traditional linear regression model.The RS data at a 20 m spatial ground-resolution came from the Airborne Visible/Infrared Imaging Spectrometer (AVRIS) in Oregon, USA.A sample of 300 pixels was taken at random from an image representing 36 km 2 .The 21st AVRIS band was selected to represent a true image.Seven other wavebands in the visible region were selected to represent RS images.The RS images obtained in these bands were treated as ancillary data.The linear correlation coefficients between the values of pixels from band 21 and the values of corresponding pixels from bands 7, 9, 10, 12, 13, 15, and 22 were 0.45, 0.55, 0.61, 0.75, 0.79, 0.88 and 0.94, respectively.Cokriging was applied to the 300 samples with each of the seven ancillary images.Without RS data, kriging would be one way to interpolate values between vegetation samples on the ground.By adding a second RS ancillary variable to a kriging analysis, it is possible to increase the accuracy of the interpolation with cokriging (Webster and Olivier 1990).The cokriging predictor at an unsampled location has the standard form: where ˆ( ) p x 0 is the prediction at location x 0, p(x j ) is the jth nearby sample value weighted by λ j , and a(x l ) is the ith nearby covariate value weighted by ω i .N 1 and N 2 are, respectively, the numbers of nearby sample values and nearby covariate values.
The accuracy in the spatial prediction of vegetation quantities was measured using the root mean square error (RMSE) between the predicted and the true image.In the case when the correlation coefficient between sample and ancillary data did not exceed 0.89, the lower RMSE was obtained using cokriging.When this coefficient exceeded 0.89, regression gave a more accurate predictor.
De Bruin (2000) presented a geostatistical method for modeling spatial uncertainty in estimates of the area extent of land-cover types.He used a sequential indicator simulation technique that enables the generation of multiple maps that take into account the available data and allows spatial patterns and uncertainties in the maps to be inferred (Deutch and Journel 1998).The area estimates were based on exhaustive but uncertain ("soft") RS data and a sample of reference ("hard") data.Using sequential indicator simulation, a set of equally probable maps was generated from which uncertainties regarding land-cover patterns can be inferred.Collocated indicator cokriging was also employed for estimation of the real extent of land-cover types.This method explicitly accounts for the spatial cross-correlation between hard and soft data using a simplified model of coregionalization.The method was illustrated using a case study from olive tree plantations in southern Spain.Demonstrated uncertainties concerned the spatial extent of a contiguous olive region and the proportion of olive vegetation within large pixel blocks.Because the image-derived olive data were not very informative, conditioning on ground-truthed data had a considerable effect on the area estimates and their uncertainties.For example, the expected area of the contiguous olive region increased from 65 ha to 217 ha when conditioning on the reference sample.
Conventional classification techniques require ground truth information, use only the spectral characteristics of classified pixels, rely on a Gaussian probability distribution for the spectral signature, and work at a pixel level without allowing classification of larger or smaller areas.To overcome these drawbacks, van der Meer (1996) proposed use of a non-parametric geostatistical technique based on indicator kriging.This method has the following advantages: 1) it can use spectral information from laboratory studies to train the classifier instead of requiring ground information, which is difficult to obtain or not always available; 2) through the kriging estimation variances, an estimate of uncertainty is derived; 3) it incorporates spatial aspects through the use of local estimation techniques; 4) it is distribution free; and 5) it may be applied to different supports if the data are corrected for support changes.Although van der Meer applied this technique to the problem of calcite-dolomite mineral mapping, this method could be promising for forestry applications assuming that reflectance spectra of the classified species type are well known and high-resolution data are used.

Conclusions
Geostatistical measures of spatial variability provide unique structural information on RS images.However, application of these measures for classification of complicated forest communities is not a trivial task.One of the most important shortcomings with the use of semivariograms is that it should be calculated only from sufficiently large and homogenous areas.Progress in RS technology allows for separate trees or their parts to be observed in the forests.The role of textural information could be therefore essential when analyzing RS images.The results obtained up to now are encouraging and the use of textural information significantly improves classification in the case of high-resolution RS data.
(McBratney and Webster 1986emivariogram is used to fit an appropriate theoretical model, such as the spherical, exponential, etc.(McBratney and Webster 1986).This model allows us to calculate semivariance values for any h that are necessary for other geostatistical calculations and analyses such as kriging or cokriging.