Accuracy assessment of a large-scale forest cover map of central Siberia from synthetic aperture radar

Can. J. Remote Sensing, Vol. 28, No. 6, pp. 719–737, 2002

Accuracy assessment of a large-scale forest cover map of central Siberia from synthetic aperture radar Heiko Balzter, Evelin Talmon, Wolfgang Wagner, David Gaveau, Stephen Plummer, Jiong Jiong Yu, Shaun Quegan, Malcolm Davidson, Thuy Le Toan, Michael Gluck, Anatoly Shvidenko, Sten Nilsson, Kevin Tansey, Adrian Luckman, and Christiane Schmullius Abstract. Russia’s boreal forests host 11% of the world’s live forest biomass. They play a critical role in Russia’s economy and in stabilizing the global climate. The boreal forests of central and western Siberia represent the largest unbroken tracts of forest in the world. The European Commission funded SIBERIA project aimed at producing a forest map covering an area of 1.2 million square kilometres. Three synthetic aperture radars (SAR) on board the European remote sensing satellites ERS-1 and ERS-2 and the Japanese Earth resources satellite JERS-1 were used to collect remote sensing data. Radar is the only sensor capable of penetrating cloud cover and imaging at night. An adaptive, model-based, contextual classification to derive ranked total growing stock volume classes suitable for large-scale mapping is described. The accuracy assessment of the Siberian forest cover map is presented. The weighted coefficient of agreement κw is calculated to quantify the agreement between the classified map and the reference data. First, the classified map is compared with Russian forest inventory data (κw = 0.72). The inherent uncertainty in the forest inventory data is simulated by allowing for fuzziness. The effect of uncertainty on the unweighted coefficient of agreement κ is stronger than that on the weighted coefficient of agreement κw. Second, the map is compared with a more reliable, independent posterior ground survey by Russian forestry experts (κw = 0.94). The follow-on project SIBERIA-II started in January 2002 and is striving to develop multisensor concepts for greenhouse gas accounting (www.siberia2.uni-jena.de). 737

Résumé. Les forêts boréales de Russie contiennent 11% de la biomasse vivante mondiale. Elles jouent un rôle primordial dans l’économie Russe et, en outre, aident à stabiliser le climat mondial. Les forêts boréales de Sibérie Centrale et Occidentale représentent la plus large étendue de forêt continue dans le monde. Le projet SIBERIA financé par la communauté Européenne a pour but de produire une carte de ces forêts sur une surface de 1.2 millions de km2. Trois systèmes de radar à synthèse d’ouverture (SAR) montés sur les satellites ERS-1, ERS-2 and JERS-1 ont fourni les données de télédétection requises pour cette étude. Le radar est le seul capteur capable de pénétrer à travers les nuages et d’imager pendant la nuit. Une procédure contextuelle de classification est développée afin d’extraire les classes de la volume de tronc requises. L’évaluation de l’exactitude de la carte de couverture forestière de Sibérie est présentée. Le coefficient pesé κw est calculé afin de quantifier la conformité entre la carte de classification et les données de référence. Tout d’abord, la carte de classification est comparée à un premier inventaire Russe des forêts (κw = 0,72). L’incertitude inhérente à cet inventaire est modélisée explicitement à l’aide d’un procédé de bruit blanc. L’effet des variations de l’incertitude est plus élevé pour le coefficient non pesé κ que pour le coefficient pesé κw. Enfin, la carte est comparée à un deuxième inventaire des forêts plus raffiné que le précédent éffectué par des experts forestiers Russes (κw = 0,94). Faisant suite à SIBERIA-I, le projet SIBERIA-II (www.siberia2.uni-jena.de), qui a débuté en janvier 2002, s’efforce à développer des concepts multi-capteurs pour l’évaluation du bilan des gaz à effet de serre.

Received 14 March 2001. Accepted 19 June 2002. H. Balzter,1 D. Gaveau, and S. Plummer. Centre for Ecology and Hydrology Monks Wood, Abbots Ripton, Huntingdon, Cambridgeshire, PE28 2LS, U.K. E. Talmon. Justus-Liebig-University, Giessen, Germany. W. Wagner. German Aerospace Research Establishment, Wessling, Germany. J.J. Yu and S. Quegan. University of Sheffield, U.K. M. Davidson and T.L Toan. Centre d’Études Spatiales de la Biosphére, Toulouse, France. M. Gluck, A. Shvidenko, and S. Nilsson. International Institute of Applied Systems Analysis, Laxenburg, Austria. K. Tansey and A. Luckman. University of Wales, Swansea, U.K. C. Schmullius. University of Jena, Germany. 1

Corresponding author (e-mail: [email protected]).

© 2002 CASI

719

Vol. 28, No. 6, December/décembre 2002

Introduction 737 The boreal forest belt of Russia is a significant carbon pool (Nilsson et al., 2000). Because of the high latitude and low temperatures, ecological and biological processes in this forest are much slower than in the tropical zone. Russia’s boreal forest is therefore currently being discussed by the carbon cycle community as a potential long-term carbon sink (Kokorin et al., 1996; Rugo and Weiss, 1996; Schulze et al., 1999). However, human-induced forest fires destroy vast areas of natural forest in Russia every year. Regrowth is slow, and the carbon release during and after a fire contributes to global warming (Kokorin et al., 1996). Because of the less developed infrastructure and long distances in Siberia, forest inventories are not carried out frequently enough to provide the information on the boreal ecosystem that would urgently be needed. Some forest inventories are more than 40 years old. To update these forest inventories optical imagery acquired by airborne and satellitebased sensors has been used operationally. In addition to field surveys and optical data, synthetic aperture radar (SAR) data are being discussed as a potential source of information about the state of the Siberian forest. In a frequently cloud covered large region like Siberia, cloud-penetrating radar sensors can provide large-area coverage at high spatial resolution. The objectives of the SIBERIA project were as follows: (i) to map a largely inaccessible area of 1.2 million square kilometres of boreal forest; (ii) to explore the potential operational capabilities of SAR for large-scale forest mapping; and (iii) to provide a spatial data product for updating total growing stock volume in Russian forest inventories for stands that have changed due to fire, logging, and regeneration. This paper focuses on the accuracy assessment of the forest cover map.

Data sets Remote sensing images Synthetic aperture radar (SAR) sensors on board three satellites were used to acquire remote sensing images. Due to a gap in the global network of satellite receiving stations, Siberia could not be covered by SAR images in the past (Schmullius and Rosenqvist, 1997). In 1997–1998, the German Aerospace Center (DLR) deployed a mobile receiving station in Ulaanbaatar, Mongolia, to record SAR signals from JERS-1, ERS-1, and ERS-2. In the SIBERIA project more than 600 JERS-1 L-band SAR images and 366 C-band ERS-1 and ERS2 SAR images were processed. The images were acquired in two campaigns, 22 September to 27 October 1997 and 10 May to 10 August 1998, the latter including JERS-1. ERS-1/2 tandem image pairs acquired during autumn 1997, ERS-2 images acquired during spring 1998, and JERS-1 images acquired during summer 1998 were mosaicked. JERS-1 carried an L-band SAR with 24 cm wavelength and horizontal transmit and receive polarization (HH), and ERS-1 and ERS-2 carry Cband SAR with 5.6 cm wavelength and vertical transmit and vertical receive polarization (VV). 720

Longer wavelengths like L-band penetrate deeper into the vegetation canopy, and scattering of the radiation originates from trunks and large branches. At shorter wavelengths like Cband, scattering occurs in the upper layers of the canopy from leaves and small branches. Slope and aspect affect radar backscatter. SAR interferometry provides information on the topography of the surface and on temporal changes in certain land surface properties. This technique can retrieve both structural information of natural targets, which in some cases can be converted to biophysical parameters, and digital elevation models (DEM), which can be used for geocoding and radiometrically correcting SAR backscatter imagery. Two SAR sensors image an area with the same sensor characteristics from similar viewing positions. For this project, ERS-1 and ERS-2 were programmed to fly a tandem mission with similar orbits separated by 1 day. The two sensors are separated by a spatial baseline. From the two signals an interferogram can be computed from which two parameters can be derived: (i) the interferometric coherence as a measure of the correlation between the two signals; and (ii) the interferometric phase, which is related to topographic height. The following two types of preprocessed ERS images were produced by the DLR: (1) Geocoded and terrain corrected (GTC) images were obtained if the data quality was high enough to derive an interferometric DEM. This DEM was then used to radiometrically and geometrically correct the image for topographic effects. The DEM production was only possible in areas of moderate topography. (2) Geocoded and ellipsoid corrected (GEC) images were produced if the coherence was too low to derive a DEM. In this case the image geometry was corrected for the Earth’s ellipsoidal shape using orbital parameters and the GTOPO30 DEM provided by the United States Geological Survey. A geometric correction but no radiometric correction for topography was carried out. All ERS images were multilooked (n = 10) and resampled to a pixel spacing of 50 m. The JERS-1 raw data were calibrated, geocoded, and geometrically corrected using the GTOPO30 DEM and fine-coregistered to the ERS frames using an automated intensity cross-correlation method by Gamma Remote Sensing. To improve the radiometric quality of all backscatter images, the data were jointly filtered using a multitemporal speckle filter, which preserves the original spatial resolution (Quegan et al., 2000a). The ERS frame and orbit numbers were used as a referencing system. Masking of high-relief areas was necessary for two reasons. First, the ERS GEC and the JERS images are not radiometrically terrain corrected. Second, the terrain-induced distortions for the different viewing geometries can make the coregistration of JERS to ERS images very inaccurate. The method developed here works as follows: (i) resample the GTOPO30 DEM to 50 × 50 m pixel spacing and coregister to © 2002 CASI

Canadian Journal of Remote Sensing / Journal canadien de télédétection Table 1. Characteristics of the test territories with forest inventory data. Location (decimal degrees) Lower left coordinates Test territory

Lat. N

Bolshemurtinskii Chunsky Ermakovsky Hrebtovskii Irbeiskii Manskii Nishne–Udinskii Primorskii Sayano–Shushensky Shestakovsky Juzhno–Baikalsky Ulkanskii Ust–Ilimsk Total

56.83 57.42 52.85 58.64 54.50 55.00 53.00 55.58 52.25 56.10 51.33 55.00 58.83

Upper right coordinates Long. E 91.83 95.17 91.48 98.37 95.25 93.00 95.83 102.09 90.50 102.94 103.08 107.75 102.67

Lat. N 57.33 58.08 53.17 59.98 55.67 55.67 55.83 55.99 53.08 56.68 51.83 55.92 59.83

Number of test areas Long. E 94.00 98.25 93.20 99.74 96.83 94.00 100.00 102.56 92.42 104.51 104.75 108.83 103.83

the respective ERS frame, (ii) calculate a geocoded incidence angle mask (GIM) based on the GTOPO30 DEM and the specific ERS acquisition geometry, (iii) calculate the standard deviation of the incidence angles for subsets of the GIM of a specific size (20 × 20 pixels gave the best results), and (iv) apply a threshold to the standard deviation to mask out hilly terrain (1.4° gave the best results). This masking method was shown to ensure the quality of the intensity images. It was validated by calculating the mean absolute difference between radiometrically terrain corrected intensity images and intensity images that were only ellipsoid corrected as a function of the threshold value and window size. Of the total imaged area of 904 000 km2, 656 000 km2 were classified and 248 000 km2 had to be masked out. Forest inventory and ground survey data Major tree species of non-mountain forests are larch (Larix dahurica and Larix sibirica) and pine (Pinus sylvestris), which cover approximately two thirds of the forested area. Larch dominates in northern regions, but is usually present in all forest formations. At the time of the image acquisitions, the leaf-shedding of non-evergreen trees was completed. Russian forest inventory (FI) data in a geographic information system (GIS) vector database provided information on total growing stock volume V (units of m3/ha, and defined as stem volume for all living species in young stands, and as all trees with diameter at breast height greater than or equal to 6 cm in old stands), land cover class, dominant tree species, topography, and other features. Total growing stock volume is the main parameter of interest. Attribute values of the FI polygons are generated by manual airphoto interpretation by the Russian forest enterprises. The FI data covered 13 test territories defined by forest enterprise boundaries (Table 1) in four administrative regions © 2002 CASI

Total 4 5 4 4 5 4 4 4 4 4 3 4 1 50

For training

4 1 2

For accuracy assessment

Training and accuracy assessment 3 5 2 1

2 4

2 2

2 4 3 1 17

1

1

7

16

Table 2. Weight matrix for the calculation of κw for the four forest classes as used in the comparison with the forest inventory (FI) data. Forest inventory data (m3/ha)

Remotely sensed data

≤ 20

20–50

50–80

>80

≤ 20 20–50 50–80 >80

1.00 0.89 0.56 0

0.89 1.00 0.89 0.56

0.56 0.89 1.00 0.89

0 0.56 0.89 1.00

of Russia (Krasnoyarsk kray and Irkutsk oblast; and small parts of the republics of Burjatija and Touva). Thirty-three test areas within these test territories were used in the classifier development. Another 12 test areas in four test territories were held back for the accuracy assessment, of which only seven were inside the imaged area and not masked for topography. These seven held-back test areas were not available to the project team during the development of the classification method. During the accuracy assessment the held-back test sites and the training test sites gave similar accuracies (see Table 4), so they were later pooled together to increase the sample size. The lack of independence of the validation dataset was compensated by an independent ground survey by the Russian foresters in 2000. The time lag between the test-site dependent updating year of the FI data and the remote sensing acquisition date introduced a problem of distinguishing classification errors from land cover change. To analyse the magnitude of this effect, an approach of spatial accuracy assessment was adopted for three test sites and is explained in the section titled Errors in the forest inventory and ground survey data. This labourintensive method was not feasible for the entire Siberian region,

721

Vol. 28, No. 6, December/décembre 2002 Table 3. Weight matrix for the calculation of κw for all six classes as used in the comparison with the ground survey (GS) data. Ground survey data Remotely sensed data

Water

Smooth open area

≤ 20 m3/ha

20–50 m3/ha

50–80 m3/ha

>80 m3/ha

Water Smooth open area ≤ 20 20–50 50–80 >80

1.00 0 0 0 0 0

0 1.00 0 0 0 0

0 0 1.00 0.96 0.84 0.64

0 0 0.96 1.00 0.96 0.84

0 0 0.84 0.96 1.00 0.96

0 0 0.64 0.84 0.96 1.00

Table 4. ERS reference frames and forest inventory test sites used for the accuracy assessment. ERS orbit

ERS frame

Type of imagea

Baseline (m)

Test site

Forest inventory updating year

κw

32543 32543 32500 32357 32657 32657 32400 32586 32500 32414 32414 32600

2439 2493 2493 2493 2493 2493 2457 2439 2403 2493 2511 2475

GEC GEC GEC GEC GEC GEC GTC GTC GTC GTC GTC GTC

230.0 244.3 247.7 273.0 169.7 169.7 219.9 187.4 224.6 227.2 233.0 180.3

Chunsky 1 Irbeiskii 2 Irbeiskii 3 Manskii Ulkanskii 1 Ulkanskii 2 Bolshemurtinskii Chunsky 2 Hrebtovskii Nishne 1 – Ukarsk Nishne 2 – Porog Primorskii

1998 1993 1996 1999 1998 1998 1998 1998 1996 1997 1997 1996

0.74 0.33 —b 0.56 0.49 0.47 0.63 0.38 0.46 0.88 0.62 0.68

Note: κw was calculated without the uncertainty model (σ ε = 0) and reflects the degree of correspondence with the forest inventory data rather than the accuracy of the map. The GIS updating year is usually 1 year after the acquisition of aerial photography. a GEC, geocoded and terrain corrected; GTC, geocoded and ellipsoid corrected. b Small n.

so a simpler approach to error modelling in the accuracy assessment was developed. After the map production, an independent posterior ground survey (GS) was carried out in June 2000 by Russian forestry experts with the aim of achieving a more reliable accuracy statistic of the map. It was motivated by the unquantified uncertainty in the FI data. Data sources used in the GS were new aerial photography, remote sensing images from other satellites (Landsat thematic mapper (TM), Satellite pour l’Observation de la Terre – VEGETATION (SPOT-VGT), advanced very high resolution radiometer (AVHRR), Earth remote sensing system (RESURS)), and data collected directly in the field. The survey was carried out with a systematic sampling scheme in seven test territories, which were selected based on the following requirements: (i) presence of major types of vegetation, landscape, land cover, and disturbances (forest fires and logging); (ii) no overlap with test areas used before by the project; (iii) recent and reliable data sources; and (iv) sample size >400 plots of 1 ha.

722

Methods The magnitude of microwave backscatter is a result of the interactions with the forest canopy (transmission, reflection, and absorption), which are influenced by the geometric and dielectric properties of the target. Three basic scattering mechanisms can be distinguished over forested areas: volume scattering, double-bounce, and surface scattering (Freeman and Durden, 1998). The surface scattering of the forest floor depends on surface roughness, periodic surface patterns, and dielectric properties (mainly water content). The volume scattering depends on the canopy thickness and density, size distribution, shape distribution, orientation distribution, and dielectric properties of the scatterers, i.e., leaves, needles, branches, and stems and on the wavelength. The L-band backscatter coefficient is generally higher for vegetated areas than for smooth dry soil and is nonlinearly related to total growing stock volume of the forest (Luckman et al., 1998; Fransson and Israelsson, 1999). Backscatter in this case is the result of the addition of canopy scattering and attenuated ground scattering. The canopy scattering and attenuation mainly come from the branches, and their sizes and orientations are of importance in determining the total backscatter. The © 2002 CASI

Canadian Journal of Remote Sensing / Journal canadien de télédétection

ground contribution is present until a high density of branches is reached that prevents radiation from reaching the ground. Overall there is an increasing trend of the L-band backscatter coefficient as a function of biomass, since the ground contribution has in general lower backscatter than the canopy for smooth soil, or very small size underlying vegetation when compared with the L-band wavelength of approximately 24 cm. The JERS-1 SAR signal is sensitive to aboveground biomass over a larger range of values than that of ERS: the L-band sensitivity is of the order of 2–3 dB for a range of biomass of about 0–70 t/ha or 0–100 m3/ha in stem volume (Le Toan et al., 1992; Luckman et al., 1997). The relationship between the JERS-1 backscatter coefficient and biomass, and also the saturation level, depends on the allometric relations between branch biomass and total biomass or stem volume. The interferometric ERS tandem coherence γ has also been observed to be associated with total growing stock (Hyyppä et al., 2000). It is defined as (Askne et al., 1997) γ = | γ | e jφ =

s1 s*2

(1)

s1 s1* s2 s*2

where s1 and s2 denote the first and second complex SAR images, here the two ERS-1 and ERS-2 images acquired with 1 day delay; j is the imaginary number; and φ is the interferometric phase. It is possible to factor the total coherence into different sources of decorrelation (Le Toan et al., 2001) such that γ = γ processing γ geometryγ volume γ temporal

(2)

(1) The decorrelation due to processing (γ processing) includes the effect of image registration and bias due to estimation of the coherence modulus. For instance, if the registration is not perfect, the coherence is reduced. The loss of coherence is accentuated on terrain with relief. (2) The decorrelation due to geometry (γ geometry) mainly concerns the interferometric baseline. The coherence decreases with the baseline. This effect is corrected during the interferometric processing by “spectral shift filtering” to normalize the coherence based on the responses of a stable surface. (3) The volume decorrelation (γ volume) characterizes the modification of the wave path inside the canopy between the two acquisitions and, to a lesser extent, the related changes in scattering mechanisms. This effect is caused by the change in the incidence angle and is small at Cband compared with the temporal decorrelation. (4) The temporal decorrelation (γ temporal) is the most important source of decorrelation in forest areas and is caused by the movements of the scatterers between two acquisitions. At C-band the scatterers are mainly needles, leaves, twigs, and small branches, which are highly © 2002 CASI

sensitive to wind effects. For different forest stands, the coherence decreases with an increase in the proportion of the scatterers in the stand. This means that, indirectly, the coherence decreases as a function of the biomass or stem volume. As the volume or biomass increase, the coherence drops and can reach the noise level. The shape of the decreasing function and the absolute levels of coherence as a function of biomass can be affected by the backscatter of the ground surface, which is a function of factors such as soil moisture and roughness. Under natural forest conditions like those in Siberia, it is expected that the intersite variation of the soil backscatter is mainly due to frame-to-frame variation of soil moisture rather than to the effect of surface roughness and underlying vegetation. Depending on the strength of the soil signal, the slope of the coherence versus biomass curve may vary. Rain and frost between the acquisitions will result in a significant loss of coherence. The resulting inversion is very site specific, as the soil conditions and the meteorological parameters must be taken into account. Another limitation for large-scale mapping is the saturation level of both the backscatter and the coherence signal. Since the coherence is more an indicator of the attenuation from the canopy than of its scattering, the saturation point does not depend on the same mechanisms as those for intensity, and saturation may occur later in favourable cases. The physics determining the coherence–biomass relationship in the case of ERS tandem data are determined by temporal and to a lesser extent volume decorrelation (Hagberg et al., 1995; Bamler and Hartl, 1998). Theoretical simulations using the coherence formulation by Bamler and Hartl (1998) showed that the effective interferometric baselines in the range of 270–170 m, like those used in the SIBERIA project, can cause decorrelation by a factor of between 0.74 and 0.84 for flat terrain and between 0.53 and 0.71 for a slope of 10°. Smaller baselines cause less baseline decorrelation. Data analysis showed that C-band coherence and L-band backscatter were related to total growing stock as quantified in the FI data, but the relation to C-band intensity was not consistent across frames (Quegan et al., 2000b). The classification algorithm is based on C-band coherence and L-band backscatter. Figure 1 shows an example of a scatterplot of C-band coherence and L-band backscatter. Water exhibits low coherence and low backscatter, smooth surfaces have high coherence and low backscatter, and forest shows a large cluster with a negative slope. The forest cluster can be partitioned to extract total growing stock volume classes. The choice of classes was made by a conference of users (Russian foresters) and the methodology team based on (i) a thorough analysis of the information content of the remote sensing data with respect to total growing stock, (ii) preliminary studies exploring 12 different classification approaches and a comparison using the weighted kappa coefficient κ w, and (iii) the user requirements for updating Russian forest inventory maps (Schmullius et al., 2001). Figures 2 and 3 show the nonlinearity of the relationship of total growing stock volume v to the two radar 723

Vol. 28, No. 6, December/décembre 2002

Figure 1. Scatterplot of coherence and JERS-1 radar backscatter. One brightness cycle is equivalent to 25% of the maximum point density. The long cluster ranges from high to low growing stock volume (left to right); the separated cluster with low coherence is open water, and the small cluster with high coherence is smooth open areas. © GIM International, reproduced from Balzter and Schmullius (2001).

Figure 2. Interferometric coherence of ERS-1 and ERS-2 in relation to total growing stock, showing the fitted exponential model. Reproduced from Wagner et al. (2000), © 2000 IEEE.

parameters. This nonlinearity and the saturation of both signals (backscatter and coherence) at a certain v made it impossible to retrieve v values higher than 80 m3/ha. Therefore all v values above that limit had to be aggregated into one large class. The saturation problem has yet to be solved by improving sensor capabilities, e.g., with polarimetric interferometry, but with the current state of knowledge, all radar-derived biomass maps have this limitation. The classes water and smooth open areas (including bogs, agriculture, and grassland) and four total 724

Figure 3. Backscatter coefficient of JERS-1 in relation to total growing stock, showing the fitted exponential model. Note that v = 0 in the FI database represents a wide range of cover types.

growing stock classes ≤20 m3/ha, 20–50 m3/ha, 50–80 m3/ha, and >80 m3/ha were defined as target classes. The uneven width of the classes is documented in Figure 7. The classes have an ecological meaning to the Russian foresters: 0– 20 m3/ha corresponds to non-forest vegetation, 20–50 m3/ha and 50–80 m3/ha to stages of forest regeneration, and >80 m3/ha to commercially important forest.

© 2002 CASI

Canadian Journal of Remote Sensing / Journal canadien de télédétection

Figure 4. Comparison of DEMs at Nishne 1 – Ukarsk. (a) ERS InSAR DEM, shaded relief. © European Commission ENV4-CT97-0743-SIBERIA, ESA 97/98, NASDA GBFM, DLR. (b) GTOPO30 DEM, shaded relief. © U.S. Geological Survey. (c) Difference between the two DEMs, enhanced. (d) Histogram of height differences.

Spatial stationarity of the estimators used for map production was a concern in the development of the classification method. The between-image variation of coherence and backscatter histograms for the forest classes was large. Adaptive models accounting for between-image variation had to be developed in response to the between-frame variation. The total growing stock classes were characterized by two exponential models with empirically derived constants (Equations (1) and (2)) (Wagner et al., 2002) from the FI data of the forest enterprises Bolshemurtinskii, Nishne–Udinskii, Chunskii, Primorskii, and Ulkanskii. The concept of these models is to find points in the histograms that correspond to the mature forest class and use these as an indicator of frame-specific shifts in the coherence and backscatter spectrum. This information is then used to modify the exponential functions according to the framespecific properties: © 2002 CASI

γ (v) = γ 75 + (0.330 + 0.581γ 75)e

−v

122.1

(3)

where γ is the ERS interferometric coherence, v is the total growing stock volume in m3/ha, and γ 75 is the smallest coherence value for which the coherence histogram of a satellite image reaches 75% of its maximum density; γ 75 was found to give a robust estimate of the left-hand ascent of the histogram of the volume class >80 m3/ha and was highly correlated with its class median at the five mentioned test sites (r2 = 0.88). The meaning of γ 75 is to quantify the location of the mature forest class in a new image without prior knowledge. Frame-specific shifts in the backscatter spectrum are included in the following equation: σ 0(v) = σ 75 − 2.46e

−v

107. 34

(4) 725

Vol. 28, No. 6, December/décembre 2002

deviations of coherence and backscatter showed little variation and were assumed constant. A Gaussian maximum likelihood classification was carried out with the model-derived means and fixed standard deviations derived from the training site analysis. Since the signatures of “water” and “smooth open areas” showed little variation between images, their means and standard deviations obtained by the analysis of training sites were used directly in the final maximum likelihood classification. The classifications were further improved by a contextual classification using the iterated contextual probability (ICP) algorithm (Baker and Balzter2). The ICP classifier algorithm is based on a Bayesian classification using adaptive a priori probabilities for polarimetric SAR data (Van Zyl and Burnette, 1992). Assuming a multivariate Gaussian distribution, the likelihood of receiving a signal vector Y with d elements given a particular class i is estimated from the spectral signatures with mean µi. Its probability density function is given by Devijver and Kittler (1982) as follows: Figure 5. Accuracy assessment of the DEM derived from interferometric SAR (InSAR DEM) using reference points of known height. The dotted line represents the regression of all points, and the dashed line the regression of all reference points ≤800 m.

where σ0 is the JERS-1 backscatter coefficient in dB, and σ75 is the maximum backscatter value for which the backscatter histogram of a satellite scene reaches 75% of its maximum density; σ75 was highly correlated with the median backscatter value of the >80 m3/ha class at the five test sites (r2 = 0.85). For the coherence model the standard error is in the order of 0.02 for the >80 m3/ha class, 0.06 for the 20–50 and 50–80 m3/ha classes, and 0.09 for the 0–20 m3/ha class. The standard error of the backscatter model is 0.22 dB for the dense forest class and between 0.49 and 0.79 dB for the lower total growing stock volume classes (Wagner et al., 2002). Figure 2 gives an example of the fitted function for γ, and Figure 3 gives an example of the fitted function for σ0 at one satellite image. The parameters γ 75 and σ75 account for the nonstationarity in terms of between-image variation of the shape of the exponential models. The coefficient of variation of γ 75 and σ75 for all classified images was 32 and 14%, respectively. The JERS-1 backscatter parameter σ75 shows less between-frame variation than the ERS coherence parameter γ 75. Mean ERS coherence and JERS-1 backscatter coefficient for the four forest classes were determined for each image by the models in Equations (1) and (2) and used as signatures in the image classification. The mean ERS coherence and mean JERS-1 backscatter varied between images; the coefficient of variation was between 17 and 25% for the coherence and between 10 and 13% for the backscatter. The standard 2

P( y| Ψ) = (2π)

−d

2 |Σ |

−1

2

 1 exp  (Y − µ i) T Σ −1(Y − µ i) 2  

(5)

where Σ–1 is the inverse of the covariance matrix; d is the number of spectral channels; y is the true backscatter vector, and Ψ is the true class to which it belongs; and T denotes the transpose. The likelihood contains the spectral information from the SAR image, including location, size, and shape of the ddimensional hyperellipsoids and the covariance between the channels. In the first iteration, uniform prior probabilities are assumed. Bayes’ theorem states that multiplication of the prior probability P(Ψ) with the likelihood P(y|Ψ) followed by normalization gives the posterior probability P(Ψ| y): P(Ψ | y) =

P(Ψ)P( y| Ψ) P( y)

(6)

where P(y) is the normalization constant. A maximum likelihood classification assigns each pixel to the class with the greatest likelihood, which is equivalent to the greatest posterior probability if the prior probability is noninformative. ICP uses a window width w, i.e., a window contains δ = w2 pixels, over which the posterior probabilities of the pixel belonging to a class k are averaged. This average is defined as the contextual prior probability in the following iteration. From the new contextual prior probability and the spectral likelihood, a new posterior probability is computed. A weighting of the contextual and the spectral information is carried out optionally by raising the contextual prior probabilities to the power β:

J.R. Baker and H. Balzter. The iterated contextual probability classifier — improvements in classification accuracy of radar images from boreal forests. In preparation.

726

© 2002 CASI

Canadian Journal of Remote Sensing / Journal canadien de télédétection

Figure 6. Autocorrelation coefficients calculated in range and azimuth directions for a mature forest stand in the Bratsk test site: (a) ERS-1 intensity correlation coefficient, 23 September 1997, mature stand 5; (b) JERS-1 intensity correlation coefficient, 4 May 1997, mature stand 5; (c) ERS tandem 80 pixel correlation coefficient in coherence, 23 and 24 September 1997, mature stand 5; and (d) ERS tandem 20 pixel correlation coefficient in coherence, 23 and 24 September 1997, mature stand 5. From SIBERIA final report, © European Commission ENV4-CT97-0743SIBERIA. β

 ∑ Pδ (Ψ | y)  i   Pi(Ψ) =  δi  nδ    

(7)

where δ i is an index denoting all pixels in a given window sized w × w, and nδ is the number of pixels in a window. The normalization constant is determined by pixel-wise summation over all products for each class: P( y) =

∑ [P( y| Ψ)P(Ψ)]

(8)

i

Adjustable parameters in ICP are the window size, the number of iterations, and the contextual weight.

Error sources The map generation is a combination of several processes with their own error sources. Geometric errors in the SAR images Some error sources involved in map production are related to the SAR imaging process: orbital stability, radiometric accuracy of the sensor, calibration (Shimada, 1999), signal-to© 2002 CASI

noise ratio, viewing geometry of the two SAR sensors for interferometry (mainly the baseline), or changing weather between acquisitions (Gens and Van Genderen, 1996; Bamler and Hartl, 1998). During image processing more errors are added: coregistration errors during resampling, bias of coherence estimation for low coherence values in the interferometric processing chain, effects of different viewing geometry of ERS-1/2 and JERS-1 on the geocoding, and topographic effects on the radiometric accuracy (Dowman, 1992). Geometric accuracy in the SAR processing chain is affected by the following: (1) Coregistration of the two single-look complex ERS SAR images for interferometric processing — Poor coregistration accuracy results in low coherence estimates. The image quality was checked by inspecting coherence histograms and scatterplots of coherence and total growing stock volume for all ERS images used for accuracy assessment. (2) FI vector registration to the ERS image — The root mean square error was usually less than 1.5 pixels. (3) Coregistration of JERS-1 to ERS images — This is discussed in more detail in the following paragraph. The 727

Vol. 28, No. 6, December/décembre 2002

exclusion of coregistration errors and ground offsets in the accuracy assessment was examined at three test sites. It increased κ w by around 0.1. (4) Terrain correction of ERS images to GEC–GTC products — As anticipated, the classification accuracy of GEC products was slightly poorer than that of GTC products. Global digital elevation data are available from the United States Geological Survey as the GTOPO30 data product. The pixel size of GTOPO30 is 30 arc seconds and thus depends on the latitude. Its height root mean square error depends on the data source that was used and is approximately 18 m for Siberia. In approximately half the study area, GTC images could be produced from the ERS tandem data, and a 50 m pixel spacing, interferometric DEM (InSAR DEM) is delivered as a by-product (Figure 4). The accuracy of the InSAR DEM was assessed by our Russian colleagues using 416 reference points of known height (Figure 5). Reference points above 800 m height were eliminated as outliers. A regression of the InSAR DEM to the reference DEM gave r2 = 0.69 (Figure 5). The fitted line has a slope of 1.04 and an offset of 53 m, which is possibly caused by different Earth ellipsoid models. The different viewing geometries of the JERS-1 and ERS satellites in conjunction with the two different DEMs generated pixel displacements in the coregistered images. Two types of coregistration were carried out: (i) JERS-1 GTOPO30 to ERS GTC InSAR DEM, and (ii) JERS-1 GTOPO30 to ERS GEC GTOPO30. A case study was conducted at test site Nishne 1 – Ukarsk to quantify the associated errors. Figure 4a and 4b show that GTOPO30 is much coarser than the InSAR DEM, which results in large differences in height for small-scale features such as little rivers and valleys. Figure 4c shows the difference between the DEMs. The InSAR DEM has smaller values than GTOPO30, which is caused by the little valleys that are only being picked up by the fine resolution InSAR DEM. The mean height of the InSAR DEM is 47 m smaller than that for GTOPO30 (Figure 4d). The ground offset or displacement this height difference causes is denoted ∆g (Schreier, 1993, p. 120): ∆g =

h tan θ

(9)

where h is the elevation height, and θ is the incidence angle. From this equation the theoretically possible ground offset between the images was estimated. Given the incidence angles in mid-range of 23° for ERS and 35° for JERS-1, the ground offset for case 1 is ∆g InSAR/ GTOPO =

728

hInSAR − htrue hGTOPO − htrue − tan 23 tan 35

(10)

where htrue is the real topographic height, hInSAR is the height of the InSAR DEM, and hGTOPO is the height of the GTOPO30 DEM. The ground offset between the two images is smaller if the errors of both DEMs have the same signs. Due to the smaller incidence angle of ERS, the image offset is more sensitive to the errors of the InSAR DEM. The ground offset for case 2 is ∆gGTOPO/GTOPO = 0.9(hGTOPO – htrue)

(11)

The ground offset between the ERS GECs and the JERS-1 images is approximately equal to the height error of the GTOPO30 DEM. For an estimated standard deviation of 25 m for the InSAR DEM and 75 m for the GTOPO30 DEM, the standard deviations of the ground offsets are 120 or 71 m, respectively; 95% of the ground offset values are in the intervals ±240 m (or ±5 pixels) for case 1 and ±142 m (or ±3 pixels) for case 2. Ground offsets are generally higher for rugged terrain than for gently undulating areas. Therefore mountainous regions were masked out for both cases of coregistration. For the remaining pixels the effect of the geometric error was reduced by a polygon erosion by two pixels at the edges of the FI forest stands to ensure that pixels from neighbouring forest stands are not displaced into other stands. Errors caused by the classification method Total growing stock is not directly measured by L-band SAR, but the backscatter signal is related to the structural and dielectric properties of the vegetation canopy. The correlation between the two ERS acquisitions decreases if a precipitation event changes the dielectric properties of the vegetation canopy and soil, or wind changes the geometry of the scatterers (movement of twigs, leaves, and needles). These effects lead to lower ERS tandem coherence. Equally, L-band backscatter increases with the moisture content in the vegetation. The classification method takes some of this variability into account by estimating scene-specific parameters; if decorrelation is strong, however, the information content of the coherence can be lost altogether. The model for retrieving total growing stock from ERS coherence and JERS-1 backscatter introduces some errors to the estimations: (i) decreasing slope of the exponential model, causing an increasing root mean square error (see also residual plots in Figure 9); (ii) unaccounted variability of γ and σ0 (Figures 2 and 3); (iii) estimation error of the class centres used in the maximum likelihood classification (see standard errors reported in the section titled Methods); and (iv) potential deviation from Gaussian distributions of γ and σ0 for each class in the maximum likelihood classification. Spatial autocorrelation Spatial autocorrelation reduces the number of statistically independent samples in a window and can significantly affect filter and analysis techniques that assume uncorrelated data. © 2002 CASI

Canadian Journal of Remote Sensing / Journal canadien de télédétection

Estimated autocorrelations in the range and azimuth directions, denoted by ρr and ρa, are plotted as functions of lag in Figure 6 for a mature forest stand in the Bratsk test site. The different plots correspond to ERS intensity, JERS intensity, and coherence images generated using both 80 and 20 pixel processing windows. The ERS backscatter data are significantly correlated, with both ρr and ρa having values around 0.5 at lag 1 and 0.2 at lag 2. The JERS backscatter data are almost uncorrelated, with correlation coefficients of less than 0.2 at all nonzero lags in both range and azimuth directions. For the coherence images, the correlation is not significant except at lag 1 if an 80 pixel window was used in the estimation, where ρr and ρa are both close to 0.5. The latter was the case in this study. Errors in the forest inventory and ground survey data The accuracy of the Russian FI data is of unknown magnitude. The estimates of total growing stock are generated by airphoto interpretation and are provided as rounded values in steps of 5 m3/ha up to 30 m3/ha and in steps of 10 m3/ha for greater values. A confidence interval as broad as ±20 m3/ha is possible for some forest stands (L. Vaschuk, personal communication). Legally, the Russian forest inventory manual requires 15% accuracy (with 95% confidence) for these growing stock estimates. This uncertainty contributes to the accuracy statistics when comparing classification results with the FI database. A second problem is temporal land cover change. Figure 7 shows a scatterplot of coherence and total growing stock for test site Nishne 1 – Ukarsk. The distribution of growing stock values shows many polygons of higher growing stock volumes, which are classified as low growing stock classes. This is

caused by the large number of polygons >80 m3/ha in the FI data and implies high expected errors of commission. Some of these polygons are classification errors, but some detect land cover change depending on the varying age of the FI data. For large time lags between inventory and remote sensing data acquisition, changes are likely to have happened on the ground, like harvesting, selective logging, forest fires, secondary regrowth, or replanting. These land cover changes are counted as misclassifications in the confusion matrix. To study the magnitude of these two processes, at the three test sites Irbeiskii, Hrebtovskii, and Lake Baikal south, polygon numbers were identified, which were thought to have changed since the airphotos to update the FI data were taken. In 100% of the cases it was confirmed by the Russian forest enterprises that management activities had been carried out or recent fires had destroyed the identified forest stands. In a classified image, isolated pixels are most likely misclassifications, but coherent patches are more likely to be clearcuts made after the FI update. Isolated pixels were eliminated during the ICP algorithm in the classification procedure if they had a radiometric signature that was not too distinct from that of their neighbours. Irregularly shaped coherent patches indicate forest fire scars. Even in the protected forests at Lake Baikal south, the identified polygons were classified as lower growing stock volume because of maintenance cutting in the stands. Due to the lack of an automated and objective method for removal of such polygons, they were not removed in the accuracy assessment. The accuracy statistics using the FI data have to be interpreted considering this effect. To undertake an independent accuracy assessment without confusion of land cover change with misclassification, the GS data were collected.

Figure 7. Simplified illustration of the errors involved in the classification: Nishne 1 – Ukarsk (orbit 32414, frame 2493). Vertical lines show class boundaries of the total growing stock classes, horizontal lines the coherence thresholds determined by the classification algorithm (source: DLR). Shaded areas show correctly classified polygons.

© 2002 CASI

729

Vol. 28, No. 6, December/décembre 2002

Figure 8. Overview of the forest cover map, a mosaic of classified SAR images of central Siberia. © European Commission ENV4-CT97-0743-SIBERIA, ESA 97/98, NASDA GBFM, DLR.

Accuracy-assessment methods Janssen and Van der Wel (1994) provide a review of accuracy assessments of land cover maps. A frequent problem in accuracy assessment is that the actual accuracy of the map cannot be determined because the reference data have associated errors of unknown magnitude. In this case only the correspondence of the map to the reference data can be quantified. A traditional nonspatial accuracy assessment only compares area estimates of land cover classes, but more recently the importance of a spatial accuracy assessment has been widely recognised (Lowell and Jaton, 2000). The coefficient κ (Cohen, 1960) was used for the accuracy assessment. It corrects for the chance agreement by estimating it from the observed marginal distributions of the confusion matrix and can be calculated from the confusion matrix as follows: κ =

730

p0 − pe 1 − pe

(12)

with the observed agreement p0 = chance agreement pe =

1 N

1 N

n

∑ p jj

and the expected

j =1

n

∑ p j• p• j , j =1

where pjC and pCj are the row and column sums, and pjj is the main diagonal of the confusion matrix. The coefficient κ only distinguishes between a correctly or incorrectly classified object by summing over the main diagonal pjj to calculate p0. For ranked classes like the total growing stock volume classes in the SIBERIA map, a modified coefficient is proposed that is weighted by the seriousness of the classification error. The weighted coefficient κ w (Cohen, 1968; Balzter et al., 2000; Gonin et al., 2000) weighs every single element pjk of the confusion matrix to determine the observed agreement p0. The weight matrix used here is based on the quadratic weighting function suggested by Gonin et al. (2000):

© 2002 CASI

Canadian Journal of Remote Sensing / Journal canadien de télédétection

Figure 9. Density plots of residuals (median v of all stands classified as class i minus standwise v from the FI database) for the four growing stock volume classes.

κw =

p0 − pe 1 − pe

(13)

1 N

n

n

∑ ∑ w jk p jk, j =1 k =1

w jk = 1 −

pe =

1 N2

n

Forest inventory data

n

∑ ∑ w jk p j• p•k,

and

j =1 k =1

( j − k) 2 (n − 1) 2

where wjk is the weight of element j,k used in the calculation of κ w. For instance, classifying a pixel of class 20–50 m3/ha in the FI data as 50–80 m3/ha is less serious than classifying it as >80 m3/ha. If many classification errors close to the main diagonal of the confusion matrix occur, κ is low but κ w indicates a better agreement because the classes are different but similar. If the weight matrix is the identity matrix, then κ w equals κ. The weighted coefficient κ w has been applied to FI data by Naesset (1996a; 1996b) and to a Landsat TM based forest classification by Foody et al. (1996). The transformation of the weight matrix used here to the matrix of disagreement used by Foody et al. is explained in detail in Bortz et al. (1990, p. 482). The first approach to assess the map accuracy uses a polygon from the Russian FI data as the basic element. As the FI data are vector data with polygon attributes and the classified map is pixel based, two methods were compared at five test sites: rasterization of polygons and tabulation of pixel counts, and aggregation of pixel classes in each polygon by calculating the median and tabulating polygon counts. The κ w values were similar, and the polygon-based method was adopted because © 2002 CASI

Results of the accuracy assessment In this section, results of the map accuracy assessment using FI and GS data are presented, and a simulation of the impact of the uncertainty in the FI data on accuracy statistics is given.

with p0 =

the map users will visually interpret the map on a forest-stand basis, not on a pixel basis. To calculate the accuracy of a classified frame, the following processing steps were carried out: (i) coregistration of FI vector database to the ERS frame using an automatic coarse registration and a manual fine registration with ground control points; (ii) masking of areas with rugged topography and polygon erosion at the edges by two pixels to reduce the impact of coregistration errors on map accuracy assessment; (iii) calculation of the median class of all pixels per polygon; and (iv) calculation of κ w (Equation (11)) with a quadratic weighting function (Table 2). Nonforest classes (water and smooth open areas) were not included in this accuracy assessment because there was an insufficient number of polygons with these classes in the FI data. The second approach compares the classified forest map with the GS data and uses clusters of four pixels (areas of 1 ha) as a basic sampling unit. The systematic sampling scheme provided samples of all six classes. The quadratic weight matrix was modified for classes “water” and “smooth open areas” to avoid any misclassification of these classes (Table 3).

The Siberian forest cover map is shown in Figure 8. It was delivered as a digital data product and as 123 hard-copy map sheets on a scale of 1 : 200 000. The hard copies serve those Russian forest enterprises that do not yet have GIS capabilities as information for sustainable forest management. Many forest enterprises do not have up-to-date spatial information on harvested areas in their forest inventories, and the secondary regrowth is not regularly monitored. A sample of 12 classified ERS reference frames was used to assess the accuracy of the forest cover map. Table 4 gives the 12 classified frames and the corresponding test territories (see also Table 1). The coefficient κ w for the FI data varies between the sites from 0.33 to 0.88. The polygon counts of the confusion matrices of all test sites were added to get a pooled confusion matrix for the overall map (Table 5, κ w = 0.72). Table 5 shows low user accuracies (high errors of commission) for the intermediate growing stock classes 20–50 and 50–80 m3/ha. This was anticipated from the high frequency of the class >80 m3/ha in the FI data in Figure 7 and the large residuals in Figures 2 and 4. The residuals of v were analysed by calculating the median v of all stands classified as class i and subtracting from it the standwise v given by the FI data. As expected from the shape of the backscatter–coherence curves in Figures 2 and 3 and the diagram in Figure 7, the density plot of the residuals spreads out as the volume class increases (Figure 9). For higher volume classes, the absolute residuals are greater. The residuals for the 731

Vol. 28, No. 6, December/décembre 2002 Table 5. Pooled confusion matrix for all test sites, comparing the map with forest inventory (FI) data. Forest inventory data Remotely sensed data

≤20 m3/ha

20–50 m3/ha

≤20 20–50 50–80 >80 Total Producer accuracy (%)

589 144 135 31 899 66

104 110 237 96 547 20

50–80 m3/ha 21 52 297 223 593 50

>80 m3/ha 136 117 1023 5327 6603 81

Total 850 423 1692 5677 8642

User accuracy (%) 69 26 18 94

Note: κ (0.43) and κw (0.72) were calculated with σ ε = 0.

Table 6. Variability of κw for all six classes across the ground survey (GS) test sites used for the accuracy assessment. Number of comparisons

κw

Plain; forest steppe zone; significantly transformed pine forests; forest inventory of 1999 based on air photography of 1998

1155

0.89

Shestakovsky

Hilly (up to 250–300 m) plain between Ilim and Kuna rivers; transformed pine and birch forests; masked areas along rivers

733

0.93

Birjusinsky

Plain; significant areas of clearcut and burns; air photography of 1998, inventory of 1999

645

0.97

Nizhne–Udinsky

Plain part; forests are transformed by fire and logging; significant areas of agricultural land and bogs

897

0.97

Igirminsky

Hilly area; basically pine and deciduous forests; inventory of 1998

900

0.96

Gremuchinsky 1

Upper terrace (plains with low hills) to the north from the Angara River; untransformed pine and deciduous forests; inconsistency generated by “patching together” of different scenes

436

0.73

Gremuchinsky 2

To the south from Gremuchinsky 1; land classes with low biomass are basically presented by bogs and burned areas

466

0.91

Forest enterprise

Description

Emel’janovsky

Table 7. Pooled confusion matrix for all test sites, comparing the map with the ground survey (GS) data. Ground survey data Remotely sensed data Water Smooth open area ≤20 20–50 50–80 >80 Total Producer accuracy (%)

Water 95

95 100

Smooth open area

≤20 m3/ha

137 19 1

20 908 76 12

157 87

1016 89

20–50 m3/ha 1 36 576 33 9 655 88

50–80 m3/ha

5 39 881 120 1045 84

>80 m3/ha

9 15 58 2182 2264 96

Total 95 158 977 707 984 2311 5232

User accuracy (%) 100 87 93 81 90 94

Note: Counts are 1 ha (4 pixels) sample plots determined by Russian forestry experts at seven test sites. κ = 0.88 and κw = 0.94.

class 0–20 m3/ha have a skewed distribution to the left because the stock volume in the forest inventory is often greater than that estimated by the class. Ground survey data The correspondence of the forest cover map to the GS data is more reliable than the FI data because (i) only reliable up-todate ground data were used, (ii) most likely these ground data did not contain any large errors, (iii) all dubious situations were discussed and checked based on initial ground data on a scale of 732

1 : 25 000, and (iv) five areas were partially inspected directly during a field trip. The weighted coefficient κ w for the seven test territories varies between 0.73 and 0.97 (Table 6). Table 7 shows the pooled confusion matrix. The pooled κ w is 0.94, and the unweighted κ is 0.88. The user accuracies for each individual class are all greater than 80%.

© 2002 CASI

Canadian Journal of Remote Sensing / Journal canadien de télédétection

Simulating the uncertainty in the forest inventory data The effect of uncertainty in the FI data on the unweighted and weighted coefficients of agreement was examined using the technique of accuracy-assessment curves (Morisette and Khorram, 2000). Uncertainty was modelled by allowing for fuzziness in the FI data. In this model the measured value V of total growing stock of a polygon in the FI data is interpreted as the centre of a fuzzy interval bounded by ±2σε containing the unknown true total growing stock η with 95% certainty: P(η ∈ [V − 2σ ε , V + 2σ ε ]) = 0.95

ε ~ N (0, σ ε )

(14)

where the white noise process ε is Gaussian distributed with zero mean and standard deviation σε, and N denotes a normal distribution. In the accuracy assessment not accounting for uncertainty in the FI data, a polygon is counted as correctly classified if the total growing stock value from the FI data falls within the class interval from the remote sensing classification. For instance, V = 15 m3/ha lies outside the class 20–50 m3/ha and is counted as a misclassification. The uncertainty model accepts a classification as correct if the interval [V ± 2σε] overlaps with the class interval. For σε = 5 and a polygon classified as 20–50 m3/ha, V = 15 m3/ha generates the interval [15 – 2·5, 15 + 2·5] = [5, 25], which overlaps with [20, 50] and is counted as correctly classified. The coefficients κ and κ w were calculated for σε of 0, 1, 5, 10, 20, 30, 40, and 50 m3/ha. The results are shown in Figure 10. For higher uncertainty both κ and κ w increase nonlinearly. The use of a weight matrix in calculating κ w makes it more tolerant to unknown uncertainty in the FI data than the unweighted κ. This tolerance of κ w against classifying a polygon as a neighbouring growing stock volume class means that in Figure 10 κ increases faster than κ w; κ and κ w tend towards similar values for high uncertainty. The slope of the increase is greatest for low uncertainty (Figure 10). For the accuracy figure from the Russian forest inventory manual, κ is approximately 0.50 and κ w is above 0.70. Accepting a higher uncertainty in the FI data (up to 20 m3/ha) gives κ values of around 0.72 and κ w values of around 0.86.

Discussion and conclusions The map presented in Figure 8 is the first radar-based forest cover map of Siberia on a large scale. During the project a number of problems with large-scale mapping using image mosaics had to be addressed. The saturation of the SAR signals (both coherence and backscatter) made it necessary to map only a small number of classes and to aggregate all growing stock values above 80 m3/ha into one large class. The between-image variation of ERS coherence and JERS-1 backscatter curves depending on growing stock volume required the development of an adaptive classification algorithm. This algorithm is not ideal in that it is based on simplified, empirically derived functions, but it served its purpose of rapidly classifying a large number of images with changing properties and without the © 2002 CASI

Figure 10. Accuracy-assessment curves for the uncertainty model of the forest inventory data. The broken line represents κw, and the solid line κ. Accepting increasing fuzziness of the ground data increases the coefficients of agreement. SD, σε in Equation (14).

required input variables to do a more sophisticated classification, e.g., using microwave radiation model inversion. Fransson and Israelsson (1999) analysed multitemporal dynamics of the ERS and JERS backscatter coefficient over Swedish boreal forest and found a difference of up to 3 dB between acquisitions for a given stem volume. The temporal variation was an additive effect constant throughout the range of stem volumes. Fransson and Israelsson suggested a calibration of the radar response using clearcut areas as reference targets. The Russian forest enterprises came to the conclusion that the forest cover map has a satisfactory quality for practical applications to their forest inventories. For the forest enterprises in Krasnoyarsk Kray and Irkutsk Oblast, the map provides an update of their FI databases for sustainable forest management. The limitation of the SAR mapping capabilities in areas of rugged topography led to about 27% masked area. From a forest manager’s perspective, however, the mountainous areas are not commercially used and have a low mapping priority. Major users of forest inventory data in Russia include the following: (i) at the local level (spatial scales from 1 : 1 000 to 1 : 50 000), managers and professionals at state forest enterprises, environment protection authorities at the district level, and private companies in the forest industry; (ii) at the regional level (scales from 1 : 50 000 to 1 : 1 000 000), regional bodies of state forest management and environmental protection, regional forest inventory and planning enterprises, regional offices of Avialesookhrana (the fire protection agency), regional governments, universities, non-governmental 733

Vol. 28, No. 6, December/décembre 2002

organizations, and large industrial forest companies; (iii) at the federal level (scales from 1 : 1 000 000 to 1 : 10 000 000), the Federal Forest Service of Russia, other federal ministries (Ministry of Natural Resources, Ministry of Extraordinary Situations, etc.), federal forest fire protection agency, universities, and non-governmental organizations; and (iv) federal agencies responsible for the compliance of Russian international commitments (e.g., resulting from the United Nations Conference on the Environment and Development (UNCED) held in Rio de Janeiro in 1992 or the Kyoto Protocol). The total carbon stored in Russian forests is estimated to be of the magnitude of 32 862 Tg (Nilsson et al., 2000) and has a global importance with respect to climate change. The forest cover map will contribute to understanding the role of the boreal forest belt in the global carbon cycle and can be used to improve the Full Carbon Account for Russia (Nilsson et al., 2000) by providing better estimates of disturbances (Gluck et al., 2000). The forest biomass of Russia has often been overestimated (Nilsson et al., 2000), and dynamic global vegetation models (DGVM) can assimilate forest cover data to make more realistic predictions. Most DGVMs do not take into account disturbances, and any data on vegetation damage compared to the potential natural vegetation can make the representation of the biosphere more realistic. The JERS-1 acquisitions for the project fill a gap in the Global Forest Monitoring Programme (GBFM) led by the Japanese National Aeronautics and Space Development Agency (NASDA) in collaboration with the National Aeronautic and Space Administration (NASA) Jet Propulsion Laboratory (JPL), the German Aerospace Center (DLR), the European Space Agency (ESA), and the European Commission Joint Research Centre. These acquisitions complement the other JERS-1 backscatter mosaics of the Eurasian and American boreal forest, Equatorial Africa, the Amazon, and southeast Asia. The SIBERIA map, however, is a multisensorderived map with more limited geographical coverage but greater thematic depth. The forest cover map also provides a basis for monitoring changes in global forest cover with future spaceborne SAR sensors like the advanced land observing satellite (ALOS) and the ESA environmental satellite Envisat. The potential for operational applications and limitations of spaceborne L- and C-band SAR interferometry for global forest mapping and monitoring was discussed. A new spaceborne Lband SAR is being planned for on board ALOS and as part of the TerraSAR mission. A new spaceborne C-band interferometric constellation is less likely, although the possibility of a small chaser satellite for Envisat or Radarsat could be investigated. The large ERS tandem data archive at ESA and the JERS-1 archive at NASDA potentially provide global coverage, and the methods developed here could be used to derive a global forest classification. As an important by-product, a DEM was generated for a large part of Siberia. The InSAR DEM can be used for terrain

734

correction and geocoding of other satellite imagery, in hydrological models, run-off prediction, and snowmelt models. The use of FI data for the accuracy assessment did not prove to be very reliable because of the unquantified uncertainty in the values of total growing stock volume. The presented uncertainty model makes it possible to work with FI data for which only very general accuracy statistics are known. Based on the model, the “true” accuracy of the map was estimated to be κ ∈ [0.52; 0.72] and κ w ∈ [0.76; 0.86] (Figure 10). An independent comparison of the map with the more reliable GS data gave a higher accuracy of κ = 0.88 and κ w = 0.94 (Table 7). Although the two approaches use a different number of classes with resulting specific weight matrices, their comparison showed a good correspondence of the FI κ w with modelled uncertainty and the GS κ w. Two conclusions are possible from these two independent accuracy-assessment exercises: (i) the uncertainty in the FI data is higher than anticipated, namely in the range of 50 m3/ha, which is where κ and κ w from the FI data reach the values from the GS data; and (ii) the accuracy-assessment procedure adopted for the FI data comparison underestimates the true accuracy because it does not account for forest fires, cutting, thinning, and regrowth in the forest stands since the forest inventory date. The user and producer accuracies are lower in Table 5 than in Table 7. Although the comparison of the map with the FI data indicated low reliability of the classes 20–50 and 50– 80 m3/ha, the GS data show very high class-specific accuracies for these classes. The high coefficients of agreement are thus not just an artefact of the high number of dense forest stands in the matrices. The follow-on project SIBERIA-II started in January 2002 and is striving to develop multisensor concepts for greenhouse gas accounting (www.siberia2.uni-jena.de).

Acknowledgements SIBERIA was partly funded by Framework 4 of the European Commission, Environment and Climate, Area 3.3: Center for Earth Observation, Theme 3: Space Techniques Applied to Environmental Monitoring (contract ENV4-CT970743-SIBERIA). SAR data are provided by the ESA third ERS Announcement of Opportunity (project AO3.120 (SIBERIA)) and NASDA JERS initiative Global Boreal Forest Mapping. The satellite data were received by a mobile receiving station of the German Remote Sensing Data Center of the DLR (DFD) at Ulaanbaatar, Mongolia. We wish to thank the other members of the SIBERIA team for their contributions to the project, namely Andrea Holz, Ursula Marschalk, Jan Vietmeier (DLR), Didier Dendal, Florence Ribbes (Centre d’Études Spatiales de la Biosphére, CESBIO), Yrjo Rauste (VTT Information Technology), Tazio Strozzi, Urs Wegmüller, Andreas Wiesmann (Gamma Remote Sensing), Norbert Etzrodt (University of Bayreuth), Hans Jonsson, Marianne Orrmalm, Roland Utsi, and Torbjørn Westin (Satellus). Wolfgang Köhler is thanked for supervising Evelin Talmon’s thesis.

© 2002 CASI

Canadian Journal of Remote Sensing / Journal canadien de télédétection

References Askne, J.I.H., Dammert, P.B.G., Ulander, L.M.H., and Smith, G. 1997. C-band repeat-pass interferometric SAR observations of the forest. IEEE Transactions on Geoscience and Remote Sensing, Vol. 35, pp. 25–35. Balzter, H., and Schmullius, C. 2001. Mapping Siberian forests. GIM International, Vol. 15, No. 1, pp. 40–43. Balzter, H., Talmon, E., Yu, J.J., Quegan, S., Gaveau, D., Gluck, M., Nilsson, S., Tansey, K., Oeskog, A., Roth, A., and Schmullius, C. 2000. Accuracy assessment of a large-scale forest map of Central Siberia derived from multi-sensor satellite-borne Synthetic Aperture Radar interferometry (InSAR). In Accuracy 2000, Proceedings of the 4th International Symposium on Spatial Accuracy Assessment in Natural Resources and Environmental Sciences, 12–14 July 2000, Amsterdam. Delft University Press, Delft, The Netherlands, pp. 31–34. Bamler, R., and Hartl, P. 1998. Synthetic Aperture Radar interferometry. Inverse Problems, Vol. 14, pp. R1–R54. Bortz, J., Lienert, G.A., and Boehnke, K. 1990. Verteilungsfreie Methoden in der Biostatistik. Springer, Berlin. Cohen, J. 1960. A coefficient of agreement for nominal scales. Educational and Psychological Measurement, Vol. 20, pp. 37–46. Cohen, J. 1968. Weighted kappa: nominal scale agreement with provision for scaled disagreement or partial credit. Psychological Bulletin, Vol. 70, pp. 213–220. Devijver, P.A., and Kittler, J. 1982. Pattern recognition. Prentice-Hall, London. Dowman, I. 1992. The geometry of SAR images for geocoding and stereo applications. International Journal of Remote Sensing, Vol. 13, pp. 1609– 1617.

retrieval of forest stand attributes. Forest Ecology and Management, Vol. 128, pp. 109–120. Janssen, L.L.F., and Van der Wel, F.J.M. 1994. Accuracy assessment of satellite-derived land cover data — a review. Photogrammetric Engineering and Remote Sensing, Vol. 60, pp. 419–426. Kokorin, A.O., Lelyakin, A.L., Nazarov, I.M., and Fillipchuk, A.N. 1996. Calculation of CO2 net sinks emissions in Russian forests and assessment of mitigation options. Environmental Management, Vol. 20, pp. S101– S109. Le Toan, T., Beaudoin, A., Riom, J., and Guyon, D. 1992. Relating forest biomass to SAR data. IEEE Transactions on Geoscience and Remote Sensing, Vol. 30, pp 403–411. Le Toan, T., Quegan, S., Davidson, M., and Martinez, J.-M. 2001. Radar remote sensing of forests: applications using satellite data. In Physical Measurements and Signatures in Remote Sensing, Proceedings of the 8th International Symposium, 8–12 January 2001, Aussois, France. pp. 33–40. Lowell, K., and Jaton, A. 2000. Spatial accuracy assessment. Land information uncertainty in natural resources. Ann Arbor Press, Chelsea, Mich. Luckman, A., Baker, J., Kuplich, T.M., Yanasse, C.D.F., and Frery, A. 1997. A study of the relationship between radar backscatter and regenerating tropical forest biomass for spaceborne SAR instruments. Remote Sensing of Environment, Vol. 60, pp. 1–13. Luckman, A., Baker, J., Honzak, M., and Lucas, R. 1998. Tropical forest biomass density estimation using JERS-1 SAR: seasonal variation, confidence limits, and application to image mosaics. Remote Sensing of Environment, Vol. 63, pp. 126–139. Morisette, J.T., and Khorram, S. 2000. Accuracy assessment curves for satellite-based change detection. Photogrammetric Engineering and Remote Sensing, Vol. 66, pp. 875–880.

Foody, G.M., Palubinskas, G., Lucas, R.M., Curran, P.J., and Honzak, M. 1996. Identifying terrestrial carbon sinks: classification of successional stages in regenerating tropical forest from Landsat TM data. Remote Sensing of Environment, Vol. 55, pp. 205–216.

Naesset, E. 1996a. Conditional tau coefficient for assessment of producer’s accuracy of classified remotely sensed data. ISPRS Journal of Photogrammetry and Remote Sensing, Vol. 51, pp. 91–98.

Fransson, J.E.S., and Israelsson, H. 1999. Estimation of stem volume in boreal forests using ERS-1 C- and JERS-1 L-band SAR data. International Journal of Remote Sensing, Vol. 20, pp. 123–137.

Naesset, E. 1996b. Use of the weighted kappa coefficient in classification error assessment of thematic maps. International Journal of Geographical Information Systems, Vol. 10, pp. 591–603.

Freeman, A., and Durden, S.L. 1998. A three-component scattering model for polarimetric SAR data. IEEE Transactions on Geoscience and Remote Sensing, Vol. 36, pp. 963–973.

Nilsson, S., Shvidenko, A., Stolbovoi, V., Gluck, M., Jonas, M., and Obersteiner, M. 2000. Full carbon account for Russia. International Institute for Applied Systems Analysis, Vienna.

Gens, R., and Van Genderen, J.L. 1996. SAR interferometry — issues, techniques, applications. International Journal of Remote Sensing, Vol. 17, pp. 1803–1835.

Quegan, S., Yu, J.J. and Le Toan, T. 2000a. Iterated multi-channel filtering of SAR images. In Proceedings of the IEEE International Geoscience and Remote Sensing Symposium, IGARSS 2000, 24–28 July 2000, Honolulu, Hawaii. Institute of Electrical and Electronics Engineers Inc. (IEEE), New York. CD-ROM.

Gluck, M., Banko, G., and Vrzal, W. 2000. Harnessing remote sensing to accomplish full carbon accounting: Workshop report. International Institute of Applied Systems Analysis, Vienna, Austria. Gonin, R., Lipsitz, S.R., Fitzmaurice, G.M., and Molenberghs, G. 2000. Regression modelling of weighted kappa by using generalized estimating equations. Journal of the Royal Statistical Society Series C: Applied Statistics, Vol. 49, pp. 1–18.

Quegan, S., Yu, J.J., Balzter, H., and Le Toan, T. 2000b. Combining unsupervised and knowledge-based methods in large-scale forest classification. In Proceedings of the International Geoscience and Remote Sensing Symposium, IGARSS 2000, 24–28 July 2000, Honolulu, Hawaii. Institute of Electrical and Electronics Engineers Inc. (IEEE), New York. CD-ROM.

Hagberg, J.O., Ulander, L.M.H., and Askne, J. 1995. Repeat-pass SAR interferometry over forested terrain. IEEE Transactions on Geoscience and Remote Sensing, Vol. 33, pp. 331–340.

Rugo, O., and Weiss, A. 1996. Preserving Russia’s carbon sink: strategies for improving carbon storage through boreal forest protection. International Journal of Environment and Pollution, Vol. 6, pp. 131–141.

Hyyppä, J., Hyyppä, H., Inkinen, M., Engdahl, M., Linko, S., and Zhu, Y.H. 2000. Accuracy comparison of various remote sensing data sources in the

Schmullius, C.C., and Rosenqvist, A. 1997. Closing the gap — a Siberian boreal forest map with ERS-1/2 and JERS-1. In Proceedings of the 3rd ERS

© 2002 CASI

735

Vol. 28, No. 6, December/décembre 2002 Symposium on Space at the Service of our Environment, 17–21 Mar. 1997, Florence, Italy. European Space Agency (ESA), ESA-SP 414, pp. 1885– 1890.

pCj

column sum of column j in the confusion matrix

pjC

row sum of row j in the confusion matrix

pjj

element in main diagonal in the confusion matrix

P(Ψ)

prior probability

P(Ψ| y)

posterior probability

P( y)

normalization constant

P( y|Ψ)

likelihood function

Schulze, E.D., Lloyd, J., Kelliher, F.M., Wirth, C., Rebmann, C., Luhker, B., Mund, M., Knohl, A., Milyukova, I.M., Schulze, W., Ziegler, W., Varlagin, A.B., Sogachev, A.F., Valentini, R., Dore, S., Grigoriev, S., Kolle, O., Panfyorov, M.I., Tchebakova, N., and Vygodskaya, N.N. 1999. Productivity of forests in the Eurosiberian boreal region and their potential to act as a carbon sink — a synthesis. Global Change Biology, Vol. 5, pp. 703–722.

r

correlation coefficient

R2

coefficient of determination

s1, s2

SAR images

T

transpose of matrix

v

total growing stock volume (m3/ha)

Shimada, M. 1999. Verification processor for SAR calibration and interferometry. Advances in Space Research, Vol. 23, pp. 1477–1486.

V

measured total growing stock volume (m3/ha)

w

window size in x and y direction

wjk

weight of element j,k used in the calculation of κw

y

true backscatter vector

Y

received signal vector, d dimensional

β

weighting factor of contextual information used by iterated contextual probability (ICP)

Ψ

(unknown) true class likelihood classification

δ

index

δi

index denoting all pixels in a given window sized w × w

ε

white noise process

φ

interferometric phase

γ

interferometric coherence

γ75

interferometric coherence value for which the coherence histogram reaches 75% of its maximum value

γgeometry, γprocessing

decorrelation as a result of geometry and processing, respectively

γtemporal, γvolume

decorrelation as a result of temporal and volume, respectively

η

(unknown) true total growing stock volume (m3/ha)

Schmullius, C., Nilsson, S., Shvidenko, A., and Vaganov, E. 2001. Russian forest inventory requirements and remote sensing parameters — operational aspects evolving from the SIBERIA Project. In Looking Down to Earth in the New Millennium, Proceedings of the ERS–ENVISAT Symposium, 16–20 October 2000, Gothenburg, Sweden. Schreier, G. 1993. Geometrical properties of SAR images. In SAR geocoding: data and systems. Edited by G. Schreier. Wichmann, Karlsruhe, Germany, pp. 103–134.

Van Zyl, J.J., and Burnette, C.F. 1992. Bayesian classification of polarimetric SAR images using adaptive a-priori probabilities. International Journal of Remote Sensing, Vol. 13, pp. 835–840. Wagner, W., Vietmeier, J., Schmullius, C., Le Toan, T., Davidson, M., Quegan, S., Yu, J.J., Luckman, A., Tansey, K., Balzter, H., and Gaveau, D. 2000. The use of coherence information derived from ERS Tandem pairs for determining forest stock volume in SIBERIA. In Proceedings of the International Geoscience and Remote Sensing Symposium, IGARSS 2000, 24–28 July 2000, Honolulu, Hawaii. CD-ROM. Wagner, W., Vietmeier, J., Schmullius, C., Tansey, K., Luckman, A., Quegan, S., Yu, J.J., Balzter, H., Gaveau, D., Le Toan, T., Davidson, M., and Gluck, M. 2002. A model-based approach for retrieving growing stock volume classes of boreal forest in SIBERIA. Remote Sensing of Environment. Submitted.

List of symbols in

the

maximum

number of classification

∆g

ground offset during coregistration (m)

h

elevation height (m)

hGTOPO, hInSAR

height of the GTOPO30 and InSAR digital elevation models, respectively

htrue

real topographic height

j

index

n

number of columns and rows in the confusion matrix

κ

coefficient of agreement

nδ

number of pixels in window

κw

weighted coefficient of agreement

N

total count in the confusion matrix

µi

signature vector for class i, d dimensional

p0

observed agreement in the confusion matrix

ρa

spatial autocorrelation in azimuth direction

pe

expected chance agreement in the confusion matrix

ρr

spatial autocorrelation in range direction

element i,k of the confusion matrix

σ0

radar backscatter coefficient (dB)

736

channels

the

d

pik

spectral

in

© 2002 CASI

Canadian Journal of Remote Sensing / Journal canadien de télédétection σ 75

radar backscatter coefficient value (dB) for which the backscatter histogram reaches 75% of its maximum value

σε

standard deviation of ε

θ

incidence angle

Σ –1

inverse of the covariance matrix in the maximum likelihood classification

© 2002 CASI

737