Large Scale Mapping of Forests and Land Cover with Synthetic Aperture Radar Data

1. Introduction

1.1. A short review on SAR for forest mapping

Forests are a key natural resource providing a range of ecosystem services like carbon sequestration, natural habitats for biodiverse fauna and flora, and providing food and fiber for human consumption. To obtain sound information for management, protection, and restoration of forests, some core information needs are: 1) mapping of forest extent, 2) identification of areas of disturbance, 3) estimation of above ground biomass or growing stock volume, and 4) estimation of stand canopy height. While the first two categories are of thematic character, hence directly detectable from remote sensing data, the latter two variables need inference from models driven with remote sensing data.

Severe storms and fire are examples of major disturbance events and remote sensing has been used operationally to identify them with low spatial resolution (> 500 m) optical imagery from sensors such as NOAA AVHRR, MODIS, ERS ATSR-1/2 and SPOT-Vegetation. While identification and rapid monitoring of disturbance events is invaluable, higher resolution sensors are needed to map areal extent of the events for resource management purposes. Typical optical remote sensors used to date for responding to such needs are carried on the Landsat-5/-8, SPOT, RapidEye, IKONOS, QuickBird, GeoEye, and WorldView satellites. With a spatial resolution ranging from 30 m to better than 0.5 m, accurate information on forest area and disturbances can be retrieved. Nonetheless, optical remote sensing is limited in areas which have significant cloud cover for long periods of the year (e.g. tropical), and in those regions where sun light is an additional limiting factor (e.g. boreal). Spaceborne Synthetic Aperture Radar (SAR) data with their cloud-penetrating and day-night measurements capability provide a key data resource, complementing optical instruments for forest monitoring.

For the purposes of estimating forest biomass and growing stock volumes at stand scale, SAR systems have been shown as a valuable mapping technique due to their sensitivity to the vertical structure of forests. First results were achieved at the beginning of the ‘90s using the three frequency polarimetric airborne AIRSAR NASA/JPL SAR system. A well-established approach for generating forest biomass maps is to relate the observed backscattering coefficient to ground measurements. Using parametric exponential regressions functions or non-parametric ensemble regression tree models, forest biomass is then estimated from the SAR intensity data. Improvements in the estimation can be achieved by combining different polarizations and/or by rationing several frequencies. In general, the radar backscatter was shown to be positively correlated with some biophysical parameters such as the aboveground biomass (AGB), tree height, tree age, diameter at breast height, and basal area. Comparison of radar data acquired at C-, L-, and P-band frequencies showed that correlation of the radar backscatter with the AGB increases with increasing radar wavelength. At these frequencies, HH- and HV-polarization provide a greater sensitivity to AGB than VV-polarization [1]. Based on hundreds of studies in different ecological regions, it has been recognized that backscatter-only approaches reach a saturation level, i.e. an increase of the radar backscatter do not correspond to an increase of the AGB. The typical saturation level observed is around 300 Mg/ha, for P-Band at 100 Mg/ha for L-Band with observations of HH and HV polarizations. To overcome the saturation problem five approaches have been pursued in the past decades:

1. Ferrazzoli et al. [3] proposes to make use of bistatic radar at L-band in a specular configuration. In order to demonstrate its feasibility, a simulation analysis has been carried out by using a microwave model of vegetated terrain. The results demonstrated that woody volume up to 900 tons/ha could be inferred, hence enabling to completely solve the saturation problem. However, this approach still remains at theoretical level, since up to date no bistatic L-band SAR systems are available.

2. By using low frequency the attenuation is significantly reduced, and the large scale structures (of the order of the wavelength) dominate the backscatter. The response from non-forested areas is therefore drastically reduced, normally much below the system noise floor. The response from forested areas, on the other hand, is dominated by the large trunk and branch structures together with coherence ground reflection interactions. Since these are generally where most of the AGB is stored, the correlation of the backscatter to AGB usually increases with decreasing frequency. CARABAS – an ultra-wideband airborne SAR system operating at VHF band (20-90 MHz) – has shown that the dynamic range of the scattering is significantly larger than at P-band (440 MHz), suggesting a greater sensitivity of the lower frequency [4].

3. Another approach takes advantage of the fact that tree height can be inferred using airborne single-pass Interferometric SAR (InSAR) dual frequency (X-and P-band) data, or alternatively, LIght Detection And Ranging (LIDAR) systems. AGB is subsequently retrieved by species using allometric equations. Moreover, by integrating into the inference function the interferometric height and the P-band backscatter at different polarizations, it has been demonstrated that the well-known saturation level could be overcome [5]. On a still experimental basis, it has been shown that canopy height (CH) can be also retrieved by using airborne single-pass L-band polarimetric InSAR technique including a forest model involving a random volume of scatterers situated over a ground scattering model [6].

4. C-band backscatter is generally deemed as useless due to the weak sensitivity with respect to biophysical properties. Moreover, the strong sensitivity to the dielectric properties of the scattering objects makes C-band backscatter unreliable for monitoring forest. For this reason, low frequencies have been preferred. Santoro et al. [7] demonstrate that accurate estimates of forest GSV can be retrieved from C-band backscatter data under the requirement that a large stack of observations is available. The estimation of the GSV is carried out by means of an algorithm, which combines hyper-temporal C-band data stacks, the inversion of a water-cloud model relating the GSV to the forest backscatter, and a multi-temporal combination of GSV estimates from each image. Traditionally, model training is based on in situ measurements for unvegetated and dense forest areas and corresponding forest backscatter measurements. The novel aspect is that these are identified by means of the MODIS Vegetation Continuous Fields (VCF) product, where the corresponding measures are computed. The hyper-temporal combination exploits the different sensitivities of the forest backscatter to GSV, which can be retrieved from the estimates of the a priori unknown model parameters. Results show that the inferred RMSE is generally below 40% for full resolution data and below 20% for aggregated versions at reduced spatial resolution. The key outcome is that the retrieved GSV is never affected by saturation, with estimates in line with in situ data up to 300 m³/ha.

5. Estimate forest biomass through regressions based on exponential function exclusively derived from single frequency single/dual polarizations, ensemble regression models at eco-regions level are developed by considering multi-sensor data and bio-and geo-physical gradient data (elevation, slope, aspect, canopy density, and land cover). Works performed in the United States and Chile have shown that ensemble regression models based on multi-temporal ALOS PALSAR-1 intensities and repeat-pass coherences at HH and HV polarization as well as Landsat data yield AGB and vegetation height accuracies with R² in the range of 0.7 to 0.85 [8].

In summary, looking at the thematic (area and disturbances) and bio-physical (AGB and CH in primis) information which can be extracted or inferred from remote sensing data, it can be stated that:

• Depending upon the eco-region, environmental conditions, and forest practices, remote sensing data should be accordingly selected – in particular considering the seasonality, i.e. vegetation phenology – and algorithms consequently adapted. For instance, due to the different practices, forest clear cuts in Amazon have a completely different response at L-band than in the tropical forest in Africa, where, typically, trunks are left on the ground for months, hence engendering a much stronger backscatter, compared to a cleared area.

• The availability of multi-/hyper-temporal SAR data considerably contributes, firstly, to enhance the data quality by significantly reducing the speckle. In this respect, it is worth mentioning that future hyper-temporal data stacks acquired from Sentinel-1A/B will play, in terms of signal processing, a relevant role for the provision of high quality data (i.e. high Equivalent Number of Looks) at highest level of detail. This will allow, at analysis level, the exploitation of a pixel based approach, simpler and less time consuming than a region based one. Secondly, the temporal component provides an additional source of information for the identification of land cover classes and the appropriate treatment of moisture related phonological influence on the backscatter variations not detectable in a single-date image. Thirdly, temporal data stacks, in particular at low wavelengths, allow the estimation of the interferometric SAR coherence, proven useful for thematic and bio-physical purposes. Also, multi-/hyper-temporal SAR data sets can additionally be used in fusion with optical data, as well as bio- and geo-physical gradient data, in order to develop specific ensemble regression models at eco-regions level for the retrieval of key biophysical parameters.

• Planned new SAR missions like ESA’s Sentinel-1 C-Band SAR and BIOMASS P-band SAR, JAXA’s ALOS-2 PALSAR, Argentina’s SAOCOM L-Band SAR, and NASA/India’s L-and S-Band NISAR will doubtlessly play a relevant role for the estimation of forest bio-physical parameters.

2. The ALOS PALSAR-1 mission

Between its launch in January 2006 and the end of the mission in April 2011, the ALOS PALSAR-1 system has acquired wall-to-wall global coverage on an annual basis, which has resulted in up to five acquisitions per year at a particular location [9]. The first-of-its-kind global observation strategy for the ALOS mission provided thus an unprecedented opportunity to take global snapshots of Earth’s natural resources at very narrow time-steps and high resolution. Figure 1 shows a pan-tropical ALOS PALSAR-1 HH/HV mosaic: around 17, 000 ALOS PALSAR-1 single-look complex data frames (coverage 70x70 km per frame) were multi-look processed to a 4-look image corresponding to 15 m pixel spacing. Subsequently, the multi-looked images were speckle filtered, radiometric calibrated, normalized and terrain geocoded using the Shuttle Radar Topography Mission (SRTM) Digital Elevation Model (DEM). Geocoded frames where finally assembled to image mosaic tiles in resolutions of 15 m, 50 m, and 100 m. A similar procedure, but in this case starting from the K&C ALOS PALSAR-1 slant range amplitude data (16 looks in azimuth and 4 looks in range, corresponding to a pixel spacing of 37.5 m in range and 50.7 m in azimuth), has been exploited by De Grandi et al. [10], where the first African ALOS PALSAR-1 HH/HV mosaic has been generated. In the image mosaic illustrated in Figure 1, the HH information channel was assigned to red, HV was assigned to green, and the ratio between the two (HH/HV) was assigned to blue. With the applied color assignment, green and yellow tones correspond to instances where both HH and HV information channels have high energy returns e.g., over forested and urban areas. Blue and magenta colors are generally found in non-forested areas, where the HH polarized energy often exhibits a higher return from the surface than the HV polarized energy.

Figure 1.

ALOS PALSAR-1 pan-tropical mosaic. Image Data © by JAXA/METI, Image Processing by WHRC/ASF/sarmap. Background image from MODIS by NASA/JPL.

L-band backscatter is correlated to increasing AGB or GSV, as, with increasing canopy density and height, the backscatter contribution from the forest floor declines and the volume scattering contribution from the canopy increases. The contribution of stem-ground interactions to the total backscatter is generally weak due to diffuse scattering at the rough forest floor and substantial attenuation of the signals in the vegetation layer unless the forest floor is flooded or the canopies are frozen. In order to capture the contribution of the environmental conditions to the measured forest backscatter, the retrieval of GSV or AGB with radar backscatter observations postulates a set of in situ measurements to tune models that relate the measured backscatter to the biophysical forest parameters to the prevailing conditions. However, extensive in situ data are often not available, either because of the vastness or remoteness of forests or because of restrictions on the use of existing measurements. Even if available, uncertainties connected to in situ measurements can be substantial. Two different approaches for model calibration that do require no or only very limited field data have been developed in recent years:

1. A number of investigators have assessed the possibility to calibrate models, relating radar observations to forest biophysical attributes, using LIDAR derived attribute estimates, which require only a limited set of in situ data for model calibration. The mapping of forest resources by means of fusion of LIDAR and SAR was tested, for instance, by Englhart et al. [11], Kellndorfer et al., [12], and Atwood et al. [14]. In [11], AGB estimates from a number of airborne LIDAR transects acquired over Kalimantan, Indonesia, were used to calibrate models, relating multi-temporal TerraSAR-X and ALOS PALSAR L-band data to AGB, and extrapolated to a 280, 000 ha area (RMSE of 79 t/ha, R² of 0.53). Kellndorfer et al. extrapolated airborne LIDAR derived estimates of CH for a 1,200 km² area in Maryland, USA, to an area of 110, 000 km² using SRTM, National Elevation Dataset (NED) and Landsat data as spatial predictor layers in an ensemble regression tree model. An RMSE of 4.4 m (Pearson correlation of 0.71) when independently validating against plot-level forest inventory data has been reported. Finally, in [14] for boreal forest, the RMSE of the AGB estimate was found to be 34.9 Mg/ha over a biomass range of 250 Mg/ha; only marginally less accurate than the 33.5 Mg/ha accuracy of the LIDAR technique.

2. Others have investigated the possibility of a fully automated algorithm that makes the retrieval based on radar backscatter data (mostly) independent of the availability of in situ data. Santoro et al. [7] presented a novel approach for the mapping of boreal GSV using multi-temporal ENVISAT ASAR ScanSAR C-band. A similar approach was presented in Cartus et al. [8] for ERS-1/2 Tandem coherence to map GSV classes in Northeast China. In both studies, the automation of the retrieval was accomplished with the aid of the MODIS VCF [13], which was used to calibrate semi-empirical models, relating the SAR/InSAR data to GSV.

In the following three sections, three case studies are presented. Particular emphasis in all these works is set 1) on the use of multi-temporal ALOS PALSAR-1 data; 2) on the data acquisition period; 3) on the wise integration with other data sets; 4) on the limitations of the ALOS PALSAR-1 data.

3. Fusion of ALOS PALSAR-1, Landsat ETM+ and ALS

Currently, small-footprint Airborne Laser Scanners (ALS) represent the most deployed type of LIDAR sensors. Numerous studies have illustrated the high performance of ALS for the estimation of forest biophysical attributes. Because of the scanning capability, ALS provide for the spatially explicit mapping of forests covered by transects of several hundred meters in width. However, wall-to-wall coverage of large forest areas with ALS is in most cases prohibitively expensive, which is why fusion with image data is required to generate wall-to-wall maps of forest attributes for larger areas. The goal of this study is to investigate robust methods for estimating CH and GSV by spatially extending ALS data using ALOS PALSAR-1 and Landsat ETM+ data, i.e. by calibrating models, relating the spaceborne data to CH and GSV, with the aid of ALS derived CH and GSV estimates. Landsat TM/ETM+ data have been considered as in several studies it was shown that a retrieval of forest biophysical parameters based on the fusion of SAR and optical data yielded higher retrieval accuracies [16, 17, 18].

The study area extended over three administrative regions: Maule, Biobio, and Araucania and covered parts of the coastal Cordillera and the Chilean Central Valley (Figure 2). The forest is dominated by even-aged plantations of Pinus radiata and to a lesser extent (<20% by area) Eucalyptus globulus. Stand-level forest inventory data for 437 stands that were collected in the timeframe of the LIDAR campaigns were provided by the ARAUCO timber company. The geospatial information company Digimapas Chile provided small-footprint airborne LIDAR data for an area of ~2.5 million ha. The data were acquired between 2006 and 2008. The used airborne platform consisted of a laser scanning system (Riegl LMS-Q560), two digital cameras (Applanix DSS 322) and navigation equipment (Applanix POS AV 401). The LIDAR data have a nominal range resolution of 2 cm and delivers an absolute vertical and horizontal accuracy of better than 15 and 25 cm, respectively. During the operation, the height and intensity of multiple discrete laser returns for each laser pulse were recorded. The laser point density on the ground varied between 1 and 3 hits/m² and the scan angles ranged up to 22.5°. Digimapas produced and delivered fully geocoded Digital Terrain Models (DTM) and Surface Models (DSM) with 1 × 1 m² pixel spacing.

Under a NASA/JAXA data agreement, the Alaska Satellite Facility (ASF) provided a multi-temporal ALOS PALSAR-1 dataset. In total, 189 Fine Beam Dual (FBD) polarization and 181 Fine Beam Single (FBS) polarization Single Look Complex images were available. The FBD data were acquired between June and December 2007 and the FBS data between January and June and November and December 2007. Ten FBD and 14 FBS images (from three different ALOS paths) covering the same area as the ALS data were used. The multi-annual ALOS PALSAR-1 acquisitions also allowed for the computation of the interferometric repeat-pass coherence, which describes the temporal stability of scattering between two images and generally decreases with increasing forest density and height. Despite the long repeat interval of 46 days – hence the increased risk of temporal decorrelation - spaceborne L-band repeat-pass coherence has shown some potential for the retrieval of forest biophysical parameters, in particular in combination with intensity measurements, when the imaging conditions were suitable. The interferometric coherence was computed for image pairs with temporal baselines of 46 or 92 days and all possible combinations of image modes (FBD-FBD, FBS-FBS, FBS-FBD). In total, coherence images for 11 acquisition date combinations were produced. The perpendicular baselines were between 40 and 900 m (Table 1). In addition, three Landsat 7 ETM+ images were obtained from the Global Land Cover Facility. Two images were acquired in December 2005 and one in January 2007 under cloud-free conditions over the study area. The L1T surface reflectance data were already calibrated and corrected for terrain as well as atmospheric effects [19].

Figure 2.

ALOS PALSAR-1 mosaic for Central Chile. RGB=HH, HV and HH/HV ratio. The white rectangles denote the areas with wall-to-wall ALS coverage and the red rectangles denote the areas used for testing the fusion of LIDAR, ALOS PALSAR-1 and Landsat ETM+.

A Canopy Height Model (CHM) with 1×1 m² pixel spacing was produced by subtracting the ALS DTM from the DSM. A suite of ALS canopy structure indices characterizing different aspects of the forests canopy structure were computed from the CHM for each stand in the inventory data. The indices comprised the percentiles of the height distribution of first returns in steps of 10%, the coefficient of variation, mean, kurtosis, skewness and several canopy density indices (i.e, the proportion of first returns from heights above different height thresholds). A Gaussian fit to the profile of first return heights was computed as an additional means of characterizing canopy vertical structure (i.e. with the number of Gaussians used to describe the profile). From the ALOS PALSAR-1 data, the average intensity and coherence were computed within each stand. From the Landsat ETM+ imagery, the average reflectance observed in bands 1 to 5 and 7 were computed.

3.1. Segmentation

Image segmentation was conducted on the ALS CHM using the multi-resolution segmentation algorithm implemented in the eCognition. Segmentation was used to delineate “stand-like” image object polygons across the entire study area (Figure 3). For the segmentation, the CHM was aggregated by means of simple block averaging from 1 to 4 m pixels as a tradeoff between preserving spatial detail and reducing the amount of data to a level that could be handled by the software in an acceptable amount of time. The segmentation parameters (i.e. scale, compactness, smoothness) were chosen so that segments had sizes comparable to the polygons in the inventory dataset (on average 8 ha) and smoothly followed stand boundaries visible in the CHM. For the segments, the same set of ALS, ALOS PALSAR-1 and Landsat ETM+ canopy structure indices were extracted as were for the inventory stands.

Figure 3.

Segmented canopy height model (left) and height profile of ALS first returns for a radiata pine stand with a GSV of 96 m³/ha and a CH of 17 m (right).

3.2. Modeling

For modeling the relationship between the suite of space- and airborne remote sensing data and the in situ measurements of CH and GSV we used randomForest [20], a data mining algorithm that has proven robust and computationally efficient and that has successfully been applied in several large-scale forest mapping efforts [17, 18]. In randomForest, a large number of regression trees are grown, each recursively partitioning the training data, considering at every node a random selection of predictors. The predictions from all regression trees are then averaged to obtain a single estimate. Each regression tree is grown using a random selection of samples. The rest of the samples, the so-called ‘out-of-bag’ cases (OOB), are estimated via the respective regression trees after training and the obtained OOB predictions for all trees are then averaged to obtain an unbiased estimate for the retrieval error.

3.3. Fusion experiments

In a two stage up-scaling approach, randomForest was used for the modeling of 1) the relationship between the ALS canopy structure indices and the in situ measurements of CH and GSV and 2) the relationship between the ALS-based estimates of GSV and CH and the ALOS PALSAR-1 / Landsat ETM+ datasets. The development of fusion models incorporating these data was performed within three 100 km² test sites. The test sites, as shown in Table 1, were selected so that 1) a wide range of stand growth stages were covered; 2) no management activities (e.g. thinning, logging, etc.) had occurred during the image acquisition timeframe; and 3) a cluster of inventory polygons (i.e. stands) was located within each site. Two of the selected test sites were located in the Cordillera along the Pacific coast and one in the Chilean Central Valley (Figure 2). In total, 105 inventory stands were located within the area of the three test sites. These 105 stands were used for validation purposes only (i.e. they were not used for the development of models, relating the ALS canopy structure indices to CH and GSV). At each of the test sites, ALS-derived estimates of CH and GSV for the segments in the ALS CHM were used as response variables in randomForest to develop new models, relating all available per-segment ALOS PALSAR-1 intensities and coherences as well as Landsat ETM+ reflectances to CH and GSV, respectively. The retrieval accuracy was assessed 1) by comparing the OOB estimates for CH and GSV to the per-segment ALS predictions of CH and GSV; and 2) by applying the developed models to the ALOS PALSAR-1 and Landsat ETM+ data extracted for the inventory polygons that were located within the area of the test sites and comparing the obtained estimates for CH and GSV to the respective in situ measurements.

TestSite	Slope(Mean/SD)	CH and GSV(mean/SD)	PALSAR-1 Data		Coherence (Bn)
Date	Mode
1	16/9°	CHFID: 29/11 m; CHALS: 25/9 m; GSVFID: 386/198 m^3/ha; GSVALS: 316/152 m^3/ha	10 Jul.; 10 Oct.; 7 Jan.; 25 Nov.; 11 Dec.	FBD; FBD; FBS; FBS; FBS	10 Jul. & 10 Oct. (503 m); 10 Oct. & 25 Nov. (304 m)
2	2.6/5°	CHFID: 19/4 m; CHALS: 17/6 m; GSVFID: 172/70 m^3/ha; GSVALS: 154/98 m^3/ha	5 Jul.; 5 Oct.; 17 Feb.; 4 Apr.; 20 Nov.	FBD; FBD; FBS; FBS; FBS	5 Jul. & 5 Oct. (597 m); 5 Oct. & 20 Nov. (214 m); 17 Feb. & 4 Apr. (567 m); 4 Apr. & 5 Jul. (320 m)
3	6.5/4°	CHFID: 21/3 m; CHALS: 16/6 m; GSVFID: 172/53 m^3/ha; GSVALS: 144/102 m^3/ha	17 Jul.; 1 Sep.; 1 Mar.; 2 Dec.	FBD; FBD; FBS; FBS	17 Jul. & 1 Sep. (42 m); 1 Sep. & 2 Dec. (887 m)

Table 1.

Properties of forest in inventory stands and segments derived from the ALS CHM (CH and GSV estimates from ALS) and the ALOS PALSAR-1 SAR/InSAR imagery available at three 100 km2 test sites.

3.4. Stand-level retrieval of canopy height and growing stock volume with ALS

Figure 4 illustrates the stand-level OOB estimates for GSV and CH when using all stand-level ALS canopy structure indices as predictors and the in situ measurements of CH and GSV obtained from 319 stands as response variables.

Figure 4.

RandomForest out-of-bag GSV and CH estimates based on stand-level ALS canopy structure indices versus in situ GSV and CH.

The retrieval accuracy is given with the coefficient of determination (R²), the root mean square error (RMSE), the relative RMSE (RMSEr) and the bias. The RMSEr represented the RMSE divided by the average GSV and CH in the in situ dataset and the bias was calculated from the difference between the average GSV and CH in the in situ dataset and the ALS predictions, respectively. In the case of the GSV retrieval, the R² was 0.81, the RMSE was 62 m³/ha and the relative RMSE was 22 % when comparing the randomForest OOB predictions to the inventory data. In the case of CH, the R² was 0.93, the RMSE was 1.7 m and the RMSEr was 7.1 %. The bias was negligible in both cases. When using independent sets of training (67 %) and testing samples (33 %), the retrieval accuracies did not differ significantly (in the range of 1 %) from the OOB results.

3.5. CH and GSV retrieval through synergy of ALS, ALOS PALSAR-1 and Landsat ETM+

Figure 5 illustrates the performance of the stand-level CH and GSV at the three 100 km² test sites when using all available spaceborne predictor layers (i.e. up to 18 layers including multi-temporal ALOS PALSAR-1 HH/HV intensities, repeat-pass coherences and Landsat ETM+reflectances) and the ALS-based per-segment estimates for CH and GSV as response variables.

Figure 5.

Estimates of GSV and CH obtained through fusion of multi-temporal dual polarization ALOS PALSAR-1 intensities and coherence and Landsat ETM+ data versus 1) LIDAR estimates of GSV and CH for segments derived from the CHM (grey dots); 2) in situ measurements of CH and GSV (black dots).

When comparing the OOB estimates of CH and GSV with the per-segment estimates from ALS, similar accuracies in terms of the R² and RMSEr were obtained at the three test sites. In the case of GSV, the R² was ~0.8 and the RMSEr was below 30% for all three test sites. In the case of CH, the R² was ~0.82–0.86 and the RMSEr was in the range of 15 to 17%. The RMSEs at test sites 2 and 3 were about 43 m³/ha (GSV) and 2.5 m (CH), respectively. At test site 1, the RMSEs were higher but since the average GSV and CH values were also higher (Table 1), the RMSEr was comparable to that obtained at the other two sites. The bias was always close to zero. The result of the independent validation using 105 inventory stands was consistent with those obtained when comparing the ALOS PALSAR-1 / Landsat ETM+ OOB estimates for CH and GSV to the respective ALS derived estimates. In the case of GSV, the R² was between 0.72 and 0.87 and the RMSEr between 15 and 25%. In the case of CH, the R² was between 0.76 and 0.86 and the RMSEr between 8 and 13%. The GSV and CH estimates for the inventory polygons generally presented a somewhat larger bias of up to 20 m³/ha and 1 m, respectively.

In order to evaluate the benefit of integrating multi-temporal ALOS PALSAR-1 FBD and FBS intensity images, repeat-pass coherence and Landsat ETM+ data, the CH retrieval has been repeated using different combinations of the spaceborne datasets. Eight different combinations of predictors were considered: (1) the best FBS intensity image (1×HH), (2) the best FBD intensity image (1×HH, 1×HV), (3) all FBD intensity images (2×HH, 2×HV), (4) all FBS/FBD intensity images (4–5×HH, 2×HV), (5) all FBD intensity (2×HH, 2×HV) and coherence images (1×HH, 1×HV), (6) all FBS/FBD intensity (4–5×HH, 2×HV) and coherence images (2–4×HH, 1×HV), (7) Landsat only, (8) Landsat and all FBS/FBD intensity (4–5×HH, 2×HV) and coherence images (2–4×HH, 1×HV). The retrieval accuracies that were achieved when using intensities from only one FBS or FBD acquisition were low with less than 50% of CH variance being explained and RMS errors in the range of 4 to 6 m at test sites 2 and 3 and 8 to 10 m at test site 1 (i.e. the test site with the highest average and maximum CH) when comparing the randomForest OOB against the ALS estimates (Figure 6).

Figure 6.

CH retrieval accuracy when using different combinations of the ALOS PALSAR-1 and Landsat ETM+ data as predictors. FBS, FBD and ETM stand for FBS/FBD intensity and the Landsat ETM+ data, respectively. FBDi and FBSi/FBDi denote the cases where intensities and coherences were used jointly. The white bars show the retrieval error when comparing the ALOS PALSAR-1 / Landsat ETM+ OOB against the ALS predictions. The grey bars refer to the comparison of the ALOS PALSAR-1/Landsat ETM+ predictions for the inventory polygons and the corresponding in situ measurements.

As was to be expected, the retrieval with the FBD images performed somewhat better than the retrieval with FBS since the FBD images included the HV intensity. The integration of multi-temporal intensity observations allowed substantial improvements of the retrieval performance. When combining all available FBD intensities, the R² and RMSE improved for about 5 to 12% and 0.25 to 0.8 m, respectively. The R² and RMSE improved further when adding the stack of FBS intensity images for 6.5 to 12% (R²) and 0.3 to 0.7 m (RMSE). The integration of the coherence images resulted as well in higher retrieval accuracies. When using all available FBS/FBD intensities and coherences, the R² and RMSE were in the range of 0.75 to 0.8 and 3 m at test sites 2 and 3 and 0.60 and 6 m at test site 1, respectively. Finally, the R² and RMSE improved significantly for 6 to 22% and 0.4 to 2 m, respectively, when adding the Landsat ETM+ data to the stack of multi-temporal intensities and coherences. When testing the retrieval with only the Landsat ETM+ data, the retrieval accuracy was roughly comparable to that achieved with the multi-temporal ALOS PALSAR-1 intensities at test sites 2 and 3. At test site 1 (the test site with the highest average and maximum CH), the Landsat ETM+ image even outperformed the ALOS PALSAR-1 imagery with only minor improvements being achieved when combining the ALOS PALSAR-1 and Landsat ETM+ imagery. The improvement of the retrieval accuracy with the successive integration of multi-temporal intensities, coherences and Landsat ETM+ was generally confirmed when comparing the randomForest predictions for the inventory polygons to the corresponding in situ measurements.

4. Mapping AGB across the Northeastern US using multi-temporal ALOS PALSAR-1 data

Based on a multi-temporal ALOS PALSAR-1 dataset comprising 655 FBD dual polarization data for the Northeastern US (Figure 7), following two topics have investigated:

• the feasibility of an automated model training and inversion approach, similar to those presented in [7, 21] for ENVISAT ASAR intensity and ERS-1/2 InSAR data in the L-band case;

• the retrieval performance at different spatial scales, considering the influence of the imaging conditions and the benefit of having multi-temporal data stacks.

Figure 7.

RGB mosaic of ALOS PALSAR-1 FBD data for the Northeastern US. The red channel shows the HH intensity, the green channel the HV intensity and the blue channel the HH/HV ratio. The black lines denote different mapping zones. The image on the right-hand side illustrates for one mapping zone the number of available FBD acquisitions.

4.1. Retrieval algorithm

Complex physically-based approaches for the modeling of backscatter as function of forest biophysical attributes have been developed that consider various aspects of the forest structure (e.g., the size and orientation of stem, branches and leaves) as well as scattering mechanisms. However, when aiming at retrieval the model formulation needs to be simple enough so that it can be inverted. A relatively simple physically-based Water-Cloud type of model [22] that has been tested extensively for retrieval with C- and L-band backscatter data models the backscatter from forest, σ⁰_for, as a sum of three contributions [23]:

${\sigma }_{for}^0=\left(1-\eta \right){\sigma }_{gr}^0+\eta {\sigma }_{gr}^0{e}^{-\alpha h}+\eta {\sigma }_{veg}^0\left(1-e^{-\alpha h}\right)$E1

The first term describes the direct backscatter from the forest floor, σ⁰_gr, through gaps in the canopy. The parameter η represents the percentage to which the ground is covered by the canopy. The second term describes the backscatter from the ground that was attenuated in the canopy. Herein, the exponential represents the two-way tree transmissivity, which is a function of the tree height, h, and the two-way signal attenuation, α. The last term describes the volume backscatter, σ⁰_veg, from an opaque canopy without gaps. The model can be re-written in the following form:

${\sigma }_{for}^0={\sigma }_{gr}^0{T}_{for}+{\sigma }_{veg}^0\left(1-T_{for}\right)$E2

where T_for represents the forest transmissivity:

$T_{for}=\left(1-\eta \right)+\eta e^{-\alpha h}$E3

The model in Eq. (1) expresses the forest backscatter as a function of height and gap fraction. According to scatterometer experiments at X- and C-bands [24], T_for can also be expressed as function of growing stock volume, GSV [m³/ha]:

$T_{for}=e^{-\beta GSV}$E4

with β being an empirical parameter. Since, in this study, AGB was the observable and GSV is commonly considered a proxy for AGB (i.e. AGB can be estimated from GSV using age-dependent expansion factors, e.g., IPCC, 2006), we replaced the volume with biomass, B [t/ha]:

$${\sigma }_{for}^0={\sigma }_{gr}^0{e}^{-\delta B}+{\sigma }_{veg}^0\left(1-e^{-\delta B}\right)$$E5

In Eq. (5), β has been replaced with δ to underline that the forest transmissivity is now expressed as function of biomass.

Two of the three unknowns in the model in Eq. (5) are related to the backscatter from open ground not covered by vegetation (σ⁰_gr) and to what is considered the backscatter from opaque forest canopies with infinite biomass (σ⁰_veg). In [7, 21], it was shown for C-band that the backscatter properties of open ground and dense forest canopies, and their temporal and spatial variations, could be identified with the aid of MODIS VCF by masking the intensity images for areas with low and high VCF canopy cover and taking the mean or median of the measured intensities in the masked areas, respectively. In the case of open ground with low canopy cover, ancillary datasets have to be used to exclude land cover classes (settlements, industrial areas, water surfaces, etc.) for which the backscatter can differ substantially from that of open ground. In the case of the intensity observed over dense forests, denoted as σ⁰_df, an additional compensation for residual backscatter contributions from the ground has to be carried out. The compensation of σ⁰_df for residual ground contributions can be accomplished with:

${\sigma }_{veg}^0\left(B_{df}\right)=\frac{{\sigma }_{df}^0-{\sigma }_{gr}^0{e}^{-\delta B_{df}}}{1-e^{-\delta B_{df}}}$E6

where B_df represents the AGB of dense forest. The estimation of σ⁰_veg requires knowledge of δ and B_df.

Once the parameters have been estimated, the model can be inverted to estimate the biomass from the SAR data. When inverting the model, some intensity measurements might in fact exceed the range of modeled intensities between σ⁰_gr and σ⁰_veg. Inversion for intensities lower than σ⁰_gr would result in negative biomass estimates, which is why a biomass of 0 t/ha has to be assigned. In the case of high intensity values exceeding σ⁰_veg, an inversion is not possible and for intensities slightly lower than σ⁰_veg, the inversion could result in biomass estimates far exceeding the range of realistic biomass values, which is why a maximum biomass level, B_max, has to be defined up to which inversion is carried out. In case of multi-temporal datasets, a weighted combination of the biomass estimates from each image covering a particular pixel location, B_i, can be computed to obtain new multi-temporal estimates, B_mt. Weights can be calculated using the difference between σ⁰_veg and σ⁰_gr, referred to hereafter as the dynamic range. The larger the dynamic range, the more weight is given to the particular biomass estimate.

4.2. AGB retrieval

Model training and inversion were carried out for each of the 1310 intensity images (655 HH and 655 HV PALSAR FBD scenes, off-nadir angle of 34.3°) available for the Northeastern United States. The Landsat-based canopy density maps and land cover maps from the National Land Cover Database [25] were used to identify sparse and dense forests in each backscatter image and to mask settlements, industrial areas, agricultural land and water surfaces, respectively. The forest transmissivity parameter δ, which describes the backscatter trend with increasing biomass, could be expected to depend on 1) the imaging conditions, and 2) the forest structural characteristics [26]. Model simulations indicated, however, that in the case of the ALOS PALSAR-1 data available for the Northeastern United States, the potential variations in the forest transmissivity parameter as function of the forest type and imaging conditions could be expected to be minor and that the use of a fixed value for δ (0.008 ha/t) for the biomass retrieval represented a justifiable compromise. For a detailed discussion of the forest transmissivity as function of AGB refer to [15]. The AGB of dense forests, B_df (Eq. 6), was determined via the plot data from the Forest Inventory and Analysis (FIA) plot network of the US Forest Service [27] with the 90^th percentile of the regional plot biomass distribution [7]. A fixed biomass offset, ΔB, was then used with respect to B_df to define the maximum retrievable biomass, B_max (B_max=B_df+ΔB), when inverting the models to estimate the biomass for all pixels in the intensity images for which the NLCD land cover map reported forest. B_max affects the retrieval results in the higher biomass ranges and therefore (primarily) derived population statistics such as county totals. That is why ΔB was adjusted so that the differences between the total biomass estimates per county and FIA county-level estimates were minimized; this was the case when setting ΔB to 30 t/ha (see below). The multi-temporal combination and mosaicing of the single image biomass maps were carried out in a single step for each of six mapping zones in the Northeastern United States. Both, HH- as well as HV-intensity based biomass estimates were considered for the multi-temporal combination. The resulting AGB map is shown in Figure 8.

Figure 8.

Forest AGB map for the Northeastern United States produced from 655 ALOS PALSAR-1 dual polarization intensity images acquired in 2007/08.

Since the ALOS PALSAR-1 data could not be compared directly to the plot data of the FIA database (note that the exact FIA plot locations are not publicly available), we compared the ALOS PALSAR-1 biomass maps to the biomass maps from the National Biomass and Carbon Dataset 2000 (NBCD) that were produced through the fusion of SRTM and NED elevation, Landsat TM and the FIA forest inventory data [17, 28]. The algorithm performance was assessed with the root mean square difference (RMSD) between the ALOS PALSAR-1 and NBCD maps calculated separately for different mapping zones (see Figure 7), which were adapted from the National Land Cover Database (NLCD) [25].

At full resolution, the RMSDs were large (>80 t/ha). This might on one side have been a consequence of noise in the ALOS PALSAR-1 data (residual speckle, small scale environmental effects) but it also has to be considered that at the level of 30 m pixels, the NBCD map contains considerable uncertainty as well; note that the NBCD map was validated at segment (i.e. hectare) level. The comparison of the maps at 30 m pixel size was therefore of limited explanatory value. Hence, the comparison at different pixel aggregation scales between 150 m and 6 km have been repeated. When aggregating the maps by means of simple block-averaging, the agreement increased substantially (Figure 9, right). The largest improvement could be observed when aggregating up to ~1 km pixel size, for which the RMSD was in a range of 20 to 25 t/ha. In Figure 9 (left), the ALOS PALSAR-1 biomass estimates for two mapping zones have been plotted against the NBCD biomass estimates for pixel sizes of 150 and 600 m. The comparison of the ALOS PALSAR-1 and NBCD biomass maps revealed a good agreement along the 1:1 line. The spread along the 1:1 line reduced substantially with increasing pixel size, which can be seen in Figure 9 (left) from the decreasing length of the vertical bars. The ALOS PALSAR-1 biomass estimates, however, tended to be lower than the NBCD estimates when approaching a biomass of 200 t/ha, indicating saturation effects in the L-band data.

Figure 9.

Left: Average (circles) and standard deviation (vertical bars) of the multi-temporal ALOS PALSAR-1 biomass estimates for intervals of NBCD biomass at zones 61 and 64. Right: Root Mean Square Difference (RMSD) between the ALOS PALSAR-1 and NBCD biomass maps as a function of the pixel aggregation scale.

4.3. Importance of multi-temporal acquisitions

To evaluate the importance of having multi-temporal stacks of L-band intensity, the retrieval performance for single images is discussed first. Figure 10 shows for the L-band intensity images covering one of the mapping zones the RMSD between the biomass estimates from single intensity images and NBCD as function of the dynamic range (at 150 m pixel size).

Figure 10.

RMSD between the biomass estimates from single HH (+) and HV (o) intensity images and NBCD as function of the dynamic range (zone 64, 150 m pixel size).

The figure clearly shows that the dynamic range can be considered one of the main factors influencing the retrieval performance. The figure also shows that the RMSDs for HV intensity images tended to be lower (40-80 t/ha) than for HH intensity images (50-100 t/ha). For a given FBD HH and HV image pair, the RMSD was always lower for the HV image (5 to 10 t/ha). The differences in the dynamic ranges were most likely a consequence of differing imaging conditions. For the images covering New York State, we compared the dynamic ranges with the weather conditions at the time of the sensor overpasses. The comparison with the weather data revealed no correlation with temperature; note that the temperature was consistently above the freezing point for all images so that no major temperature related fluctuations of the dielectric properties of the trees were to be expected. Weak negative correlations were observed when relating the dynamic range to the total amount of rain in the days prior to the sensor overpasses. In both polarizations, there was a trend towards lower dynamic ranges with increasing amounts of precipitation (i.e. with increasing wetness of the soils and vegetation). The Pearson correlation coefficients were between -0.3 and -0.5 depending on which timeframe prior to the sensor overpasses was considered. This result is consistent with previous finding [29].

Figure 11 demonstrates the benefit of having multi-temporal data for an area where five FBD images were available. The dashed line shows the RMSD for each intensity image (HV: 55-59 t/ha, HH: 58-71 t/ha, 150 m pixel size), the solid line shows the change in RMSD when successively integrating the particular images into the multi-temporal retrieval.

Figure 11.

Effect of the multi-temporal combination on the agreement of the ALOS and NBCD biomass maps (at 150 m pixel size). The circles connected by the dashed line denote the RMSD for each image and those connected by the solid line show how the RMSD changed when successively integrating the single image estimates into the multi-temporal retrieval.

The multi-temporal combination resulted in a clear improvement of the RMSD when combining the available 5 HV images for about 10 t/ha (compared to the best HV image). When, in addition, integrating the corresponding HH images, only slight additional improvements of the RMSD for about 3 t/ha were achieved. These results confirmed that the multi-temporal combination allowed for significant improvements of the biomass estimates. However, it has to be noted that the multi-temporal coverage acquired by ALOS PALSAR-1 was not consistent across larger areas as regionally between one and five images were acquired per year. As a result of the varying multi-temporal coverage (see Figure 7), the accuracy of the map in Figure 8 was not consistent across the entire study area. The comparison with NBCD confirmed that locally, the performance of the retrieval depended strongly on the multi-temporal coverage (i.e. the number of images) as well as the weather conditions under which the particular set of images have been acquired. The results therefore stress the need for consistent multi-temporal acquisition strategies for upcoming spaceborne L-band missions.

4.4. Accuracy at county scale

The comparison of the ALOS PALSAR-1 biomass maps with NBCD indicated that, at least at aggregated scales, the spatial distribution of biomass could be depicted with the retrieval approach presented. To further assess the performance of the retrieval algorithm, ALOS PALSAR-1 biomass maps have been compared to FIA county-level total AGB statistics for 143 counties in Maine, New Hampshire, Vermont, Massachusetts, Rhode Island, Connecticut, New York and New Jersey. The county statistics were obtained via the EVALIDator online inventory tool of the US Forest Service. For the comparison of the ALOS PALSAR-1 and FIA estimates, the ALOS PALSAR-1 per-pixel (30x30 m²) biomass estimates were summed per county. The comparison of the ALOS PALSAR-1 and FIA county-level estimates of total AGB resulted in a coefficient of determination (R²) of 0.98 and a root mean square error (RMSE) of 2.75 10⁶ t. When calculating the average biomass per county (i.e. the total biomass divided by the county size in hectares), the RMSE was 12.9 t/ha and the R² was 0.86 (Figure 12).

Figure 12.

Average AGB according to the ALOS PALSAR-1 biomass map and FIA biomass statistics for 143 counties in the Northeastern United States.

The FIA county statistics included information about the sampling error, which could be used to approximate the confidence intervals of the FIA estimates [27]. The sampling error was between 2 % and 110 % for the largest counties (with the largest number of sample plots) and for the smallest counties, respectively; the size of the counties ranged from 64 to 17, 686 km². The ALOS PALSAR-1 total AGB estimates were within the 95 % confidence intervals of the FIA estimates for most (92 %) of the 143 counties.

5. Synergetic use of multi-temporal ALOS PALSAR-1 and ENVISAT ASAR data for forest and agricultural mapping at national scale in Africa

The forest cover in Malawi is approximately 32, 000 sqkm corresponding to 34% of country surface [30]. Natural forests represent the remainder of the Miombo (Swahili word for Brachystegia) forests that once covered almost the whole country. The area of natural forests (primary forest 29%, other naturally regenerated forest 60%) over the years has remained unchanged – annual change rate between 2005 and 2010 is less than 1% – with the exception of forest reserves that have continued to grow in number. Characteristically, the trees shed their leaves for a short period in the dry season (June to October) to reduce water loss, and produce a flush of new leaves just before the onset of the rainy season (November-December to March-April). Planted forests consist of softwood (mainly Pinus patula) and hardwood species (mainly Eucalyptus).

Thanks to the availability of a multi-year ALOS PALSAR-1 FBD data set acquired in a systematic way between 2006 and 2011 [9] and a seasonal ENVISAT ASAR Alternating Polarization (AP) data stack planned and regularly acquired during the wet season, the possibility is explored to generate countrywide a forest map and a cultivated area map during the crop wet season. In synthesis, the method is based on the synergetic use of multi-temporal data by considering the different data characteristics and given acquisition modes, the acquisition periods, and the vegetation phenology during the acquisitions. Moreover, since the ultimate goal is to provide this type of information over large areas, i.e. at least national scale, and to upscale the proposed method to other regions [31], the processing chain has been developed in such a way that the products are generated automatically. Finally, it is worth mentioning (cf. Fusion of interferometric SAR and photogrammetric Elevation Data) that using the same ALOS PALSAR-1 repeat-pass InSAR data, a DEM with higher quality than the SRTM one can be provided in those nearly equatorial regions characterizes by non-dense forest condition.

5.1. Method

The data processing flow can be divided into two distinct steps. The first one converts the multi-temporal SLC data into terrain geocoded backscattering coefficient (σ°). In addition, for the 46-days ALOS-PALSAR-1 and 1-day Cosmo-SkyMed Stripmap interferometric data, coherence (γ) is computed and terrain geocoded. Noteworthy, but not discussed here, are the ionospheric effects in the equatorial region observed in the L-band data [35]. In this case, less than 10% of the ALOS PALSAR-1 images have been omitted from the processing due to severe (several dB) and non-systematic striping along the azimuth. However, in other cases [31] around one fourth of the ALOS PALSAR-1 data could not be used. In the second step, the forest and the cultivated area map are generated.

1. Sigma nought (σ°) is derived as follows:

   • Co-registration – Images acquired with the same observation geometry and mode are co-registered in slant range geometry. This step is mandatory to allow time-series speckle filtering.

   • Time-series speckle filtering – Within the multi-temporal filtering an optimum weighting filter is introduced to balance differences in reflectivity between images at different times [32]: this allows enhancing significantly the radiometric resolution and preserve the spatial one. Multi-temporal filtering is based on the assumption that the same resolution element on the ground is illuminated by the radar beam in the same way, and corresponds to the same slant range coordinates in all images of the time series. The reflectivity can change from one time to the next due to a change in the dielectric and geometrical properties of the elementary scatters, but should not change due to a different position of the resolution element with respect to the radar.

   • Terrain geocoding, radiometric calibration – A backward solution by considering a DEM is used to convert the positions of the backscatter elements into slant range image coordinates. The three dimensional object coordinates given in a cartographic reference system into the two-dimensional row and column coordinates of the slant range image are transformed using a range-Doppler approach. Radiometric calibration is performed at the same time by means of the radar equation, where scattering area [33], antenna gain patterns and range spread loss are considered.

   • Radiometric normalization – In order to compensate for the range dependency and topographic effects, the backscattering coefficient is normalized according to a modified cosine law by considering the relationship between the local incidence angle and the backscattering coefficient of forest [14] and agriculture, respectively.

   • Anisotropic Non-Linear Diffusion Filtering – This filter significantly smoothes homogeneous targets, whilst also enhancing the difference between neighbouring areas. The filter uses the diffusion equation, where the diffusion coefficient, instead of being a constant scalar, is a function of image position and assumes a tensor value [34]. In this way, it is locally adapted to be anisotropic close to linear structures such as edges or lines.

   • Removal of atmospheric attenuation – Although microwave signals have the ability to penetrate through the clouds, in case of severe storms, as it often occurs in Africa during the raining season period, it may occur at C-band that locally the backscattering coefficient is attenuated by water vapor in the range of several dB. The temporal signature of the backscatter coefficient can be affected in two ways: i) the thick layer of water vapor generates a strong decrease of backscattering coefficient, followed by a strong increase; ii) the strong rainfall generates a strong increase of the backscattering coefficient, followed by a strong decrease. These effects are corrected in the processing chain by analyzing the temporal signature: anomalous peaks are identified and the backscattering coefficient values are corrected by means of an interpolator. The correct application of this process relies strongly on a priori knowledge of the land cover type and the weather conditions when the image was acquired.

1. Forest area – Given the regular country-wide ALOS PALSAR-1 FBD acquisitions between 2006 and 2011, the high acquisition rate between June and October (dry season), and that forest extent in Malawi is relatively constant (annual change rate between 2005 and 2010 is less than 1%), a pseudo-seasonal time-series with a relatively high temporal occurrence by combining all years of ALOS PALSAR-1 FBD observations into one year has been developed. The seasonal HH-HV forest signatures – except some rare outliers representing clear cuts – are typically characterized by a constant high σ°, while other cover types show noticeable lower (bare soil, water, dry or low vegetation) or higher (settlements) values. It has to be pointed out that this net discrimination between the different land covers is given by the combination of two factors. First, the selected dry season period, which, due to the very limited vegetation development, considerably reduces the confusion between low forest biomass and others vegetation types, agriculture (i.e. maize) in particular. Second, the long wavelength, on one hand tends to smooth out the radar backscatter from land covers characterized by limited surface roughness, on the other hand it guarantees a relatively high backscatter of those land covers, like forest, where double bounce and volume scattering prevail. Note that this separability is not obtainable at shorter wavelengths, due to similarities in the volume scattering contributions.

2. Cultivated area – Pre-requisite for the generation of this product is to obtain a seasonal data set as far as possible regularly acquired along the whole crop season and according to the crop practices. This allows i) to reduce the confusion between cropped areas and the surrounding vegetated (non crop) areas; and ii) to monitor the crop development [36, 37]. The specific sensitivity of active microwave short wavelength sensors (C- and X-band) to soil properties, such as roughness and moisture content, enables the possibility to detect these changes already at the earliest stage during the field preparation, i.e. ploughing (high backscatter) and harrowing (low backscatter). During the second phase (namely from flowering to plant drying stage), the dielectric and structure properties of the plant are the key factors determining the high reflectivity at C-band HH polarization. Finally, the lower reflectivity during the plant drying is caused by the loss of plant moisture. As proposed in [36, 37], an efficient way to quantitatively describe the σ^o temporal crop signature is to derive appropriate temporal features, i.e. the relative minimum and maximum including corresponding dates; their difference; the minimum and maximum ratio between two subsequent acquisitions. These features are used to generate the product according to:

   1. The crop start of season, which is identified when there is a relative minimum followed by a maximum increment between two subsequent acquisitions;

   2. The crop peak of season, which is identified when there is a relative maximum followed by a minimum increment between two subsequent acquisitions;

   3. The pixel, which is classified as crop if:

      • conditions 1 and 2 are satisfied;

      • the range between relative minimum and maximum attains a minimum value;

      • the temporal duration between 1 and 2 is within a given duration.

5.2. Forest and cultivated area products

To cover the entire country, 65 ALOS PALSAR-1 FBD/FBS frames distributed over 5 adjacent tracks are necessary. Three to four coverages per dry season per year during four years have been used, resulting in the around 900 SLC scenes. A coherence-intensity mosaic using 46-days interferometric ALOS PALSAR-1 FBS data acquired during the wet season (January-February 2008) has been additionally generated. A total of 225 ENVISAT ASAR AP images acquired five times from October 2007 to April 2008 are used to monitor the crop growth along the whole wet season. The step to transform the SLC data into geo-referenced σ° and γ products is doubtless the most time consuming one. For this reason, the data processing is performed by means of a PC based cluster solution. Noteworthy, the algorithms have been written in a such way to fully exploit the characteristics of the processors. This setting allows to carry out the processing without any intervention of the operator and to process the data sets within few days.

Figure 13 and 14 illustrate four multi-temporal mosaics: the multi-year ALOS PALSAR-1 FBD mosaic acquired during the dry season (Figure 13, left); the seasonal (October to April) ENVISAT ASAR HH mosaic (Figure 13 right); an ENVISAT ASAR HH mosaic (October and January, the months showing the most significant radiometric changes) combined with the ALOS PALSAR-1 HV July one (Figure 14 left); the ALOS PALSAR-1 HH coherence-intensity mosaic (Figure 14 centre).

Figure 13.

(left) Multi-year ALOS PALSAR-1 FBD mosaic, 15m acquired during the dry season (mean HH=red, mean HV=green, mean HH / mean HV=blue); (right) Seasonal ENVISAT ASAR HH mosaic, 15m (ASAR HH October=red, ASAR HH December=green, ASAR HH January=blue).

Figure 14.

(left) ENVISAT ASAR HH-ALOS PALSAR-1 HV mosaic, 15m (ASAR HH December=red, PALSAR-1 HV July=green, ASAR HH January=blue); (centre) ALOS PALSAR-1 HH (FBS) coherence-intensity mosaic, 10m (coherence January-February=red, mean intensity=green, intensity difference=blue); (top right) detail of figure left; (bottom right) detail of figure centre.

These examples clearly show that depending upon the selected sensor, the acquisition mode and time, and a smart data integration, different types of information (i.e. products) can be derived. It is worth mentioning that the purpose of data synergy is not exclusively to obtain higher accuracies or more information, but it is also intended to simplify and automatize the products generation. Conditio sine qua non is, on one hand, to recognize the sensor capabilities and limitations – including the processing techniques – on the other hand, to understand the object, its characteristics, and the environmental surrounding conditions.

The multi-year ALOS PALSAR-1 FBD mosaic (Figure 13 left) acquired during the dry season undoubtedly shows a clear distinction between forest (green colour) and other cover types (blue tonalities). Note that this net separation is the outcome of i) the excellent quality of the multi-temporal speckle filtering, and ii) the averaging (omitting some rare outliers) of the HH intensities, and the HV respectively over the dry season period and years (total of 14 images). Both operations contribute to considerably improve the signal-to-noise ratio and to enhance the level of detail. Concerning the temporal averaging, it has to be pointed out that this procedure is more than reasonable, because the SAR data have been exclusively selected in the dry period, where the vegetation phenology is stable, in terms of roughness and dielectric properties, and, the forest extent variations are negligible. In order to demonstrate the usefulness of this approach, for a selected area, a 1-day interferometric Cosmo-SkyMed Stripmap image pair (3m resolution) acquired in September has been acquired and processed. Figure 15 shows, for comparison purposes, the single-date and multi-year ALOS PALSAR-1 products, and the 1-day interferometric Cosmo-Skymed Stripmap pair (3m resolution). The high level of detail of the multi-year ALOS PALSAR-1 product (bottom left) is particularly appreciable by comparing the speckle effect in the forest in the two ALOS-PALSAR-1 products, and the better feature delineations (for instance in the riparian forest) in the multi-year one. Furthermore, by visually inspecting the single trees in the multi-year ALOS PALSAR-1 product and in the Cosmo-SkyMed one, it is noticeable that almost all single trees identifiable in the 3m image are easily recognisable in the 15m multi-year scene, hardly in 15m single-date one. With respect to the use of X-band data for forest applications, it should be shortly remarked that, in general, the 1-day Cosmo-SkyMed coherence or the bistatic TerraSAR-X-Tandem one are essential, because the limited σ^o dynamic range at this frequency significantly reduces the discrimination capabilities in vegetated areas, if exclusively intensity is used. For details refer to Holecz et al. [38].

Figure 15.

(top left) Single-date ALOS PALSAR-1 FBD mosaic, 15m (mean HH=red, mean HV=green, mean HH / mean HV=blue); (bottom left) Multi-year ALOS PALSAR-1 FBD mosaic, 15m; (bottom right) Cosmo-SkyMed Stripmap coherence-intensity, 3m (1 day coherence=red, mean intensity=green, intensity difference=blue).

As expected, due to the shorter wavelength, forest is quite less distinguishable in the seasonal ENVISAT ASAR mosaic (Figure 13 right) confirming the conclusions presented in the comparative study carried out by Mitchell et al. [39]. The predominant colours in this image are the cyan (southward) and blue-violet (northward), corresponding to the crop growth after the start of the rain (around November). This is reflected in the December (green) and January (blue) acquisition, where the C-band backscattering coefficient almost attains its highest values due to the rapid crop growth (i.e. increase in the surface and volume scattering) and the high crop moisture content (increase in the dielectric constant). It is interesting to note the correspondence between the blue areas in the ALOS PALSAR-1 FBD mosaic (bare soil during the dry season) becoming brightly coloured in the ENVISAT ASAR mosaic (cultivated area during the wet season). This is becoming obvious in Figure 14 left (and detail, top right), where the ALOS PALSAR-1 HV July mosaic (green) has been merged with the ENVISAT ASAR HH mosaic of December (red) and January (blue). This colour composite demonstrates, in a qualitative but evident way, that data synergy is undoubtedly beneficial, if exploited in a wise manner.

Coherence is unquestionably a complementary and valuable source of information. However, in practice, useful interferometric correlation is often not easy to obtain due to unfavourable baseline conditions and unsystematic atmospheric effects. Moreover, the temporal decorrelation sometimes causes interpretation uncertainties. Figure 14 centre (and detail, bottom right) shows a 46-days repeat-pass ALOS PALSAR-1 HH coherence-intensity mosaic generated from an FBS image pair acquired in January and February 2008. As expected, forest has a low correlation and an average high intensity (green), meaning that the forest didn’t vary in this gap of time; blue (large intensity difference), particularly visible in the detail bottom right, represents the growing fields; the red tonalities correspond to those areas where the cover changes were minimal, i.e. meaning primarily rough bare soil areas. Nonetheless, it should be considered that tiny crops at L-band are almost transparent, hence resulting into a relatively medium coherence, as the rough bare soil. It turns out that the cultivated area is underestimated in favour of bare soil. In synthesis, a forest map could be generated even if, due to the long repeat-pass and the different baselines, the uncertainties can be relevant; these ambiguities are significantly higher for the cultivated area also considering that for this product long acquisition time intervals are unsuitable [37]. A final example on coherence is illustrated in Figure 16, which it has been obtained from a July-August FBD image pair.

Figure 16.

(left) ALOS PALSAR-1 46-days HH July-August coherence; (centre) HV July-August coherence; (right) multi-year ALOS PALSAR-1 FBD image, 15m resolution (mean HH=red, mean HV=green, mean HH / mean HV=blue).

Observing the coherence, the Miombo forest is not distinguishable in the HH polarization (average γ is 0.6), hardly detectable in the HV (average γ is 0.5). In the multi-year ALOS PALSAR-1 FBD intensity colour composite, as extensively presented and discussed above, forest is clearly separable from the surrounding land covers. It is anticipated that the perpendicular baseline of the interferometric pair is 280m, which for this frequency, is thereby appropriate for thematic analysis. The Miombo forest – bare in the dry season – on average has a tree height significantly less than 10m and a diameter breast height ranging from 10 to 20 cm corresponding to a relatively low biomass (often considerably less than 100 tons/ha). This means that the main scattering contribution is the volume, primarily induced by the tree branches; hence the HH radar backscatter is markedly attenuated. This is reflected in both, HH coherence (Figure 16 left) and HH intensity (Figure 16 right). At HV intensity, forest is well recognizable in green (Figure 16 right). This distinction is mainly given by the temporal averaging over the dry seasons and four years, which strongly enhances the separability of the various land types. On the opposite, for the HV coherence (Figure 16 centre), the discrimination between forest and the surrounding area (both having a relatively high coherence) is very limited: this is due, on one hand, to 46-days temporal decorrelation, on the other hand, to the impossibility to perform a temporal averaging. In summary, coherence is doubtless a valuable source of information, however, it should be used with care.

Based on the above considerations and evaluations, the most suitable solution is to use the multi-year ALOS PALSAR FBD data set for the forest area product, and the seasonal ENVISAT ASAR one for the cultivated area. Furthermore, in order to understand the contribution of the seasonal ENVISAT ASAR for the forest area (Figure 17 left), the obtained cultivated area (Figure 17 right) is merged in IF condition with the forest area one.

The two maps are performed using a hierarchical prior knowledge-based classifier. For the forest area the input data were the mean HH intensity and the corresponding mean HV intensity of the multi-year ALOS PALSAR FBD data set. Concerning the cultivated area, the temporal features of the seasonal ENVISAT ASAR data set have been used.

5.3. Accuracy at national scale

Validation involves the collection of ground reference data for the validation of remote sensing based products. Usually it is carried out by sampling units, i.e. points unambiguously identified by co-ordinates. As shown in Figure 18, systematic grids are used with randomly selected starting corner co-ordinates in order to ensure a representative and spatially well-distributed sample.

Figure 17.

(left) Forest area generated from multi-year ALOS PALSAR-1 FBD data; (right) Cultivated area generated from seasonal ENVISAT ASAR data.

From an operational perspective, the grouping of sample points in clusters is encouraged. Although this approach may introduce some degree of statistical bias, by significantly reducing travelling time between sample locations, it leaves more resources available for data collection in the field. The relevant parameters of this systematic cluster approach are: 1) distance between clusters; 2) number of points per cluster; 3) distance of points within cluster. The values can be fixed considering several criteria and constraints: 1) available budget; 2) size and shape of the area of interest; 3) resolution of the remote sensing images; 4) logistics; 5) average dimension of the area to be classified. In this case, it has been opted for 1km distance between the clusters, 16 points per cluster, and 250m distance between the points within the cluster. A total of 868 valid points have been collected.

Table 2 shows the obtained confusion matrices for the forest product exclusively based on the multi-year ALOS PALSAR-1 FBD (top) and for the forest product generated by combining the forest and cultivated area product (bottom). In general, the obtained accuracies are high. The main difference between the two tables is in the reduction of the omission errors of the classes sugar cane and crop, and urban. Sugar cane is a tall crop with a long crop season (typically one year). This fully explains the large omission error. By merging the cultivated area map, this error could be reduced by 30% only, because the cultivated area product, as defined here, exclusively considers the crop planted at the start of the rainy season. The same explanation is valid for the class crop. In essence, in order to almost completely remove this omission error, an annual (and not just seasonal) ENVISAT ASAR monitoring should be carried out. Concerning the class urban (mainly small cabins in the rural areas), the omission error has been reduced by 40%, leading to a 10% error. In this specific case, the radar backscatter is often random, therefore the combination of the two frequencies strongly contributes to better detect this land cover type.

Figure 18.

Validation scheme.

With respect to the cultivated area, the obtained overall accuracy is 80%. For details refer to Holecz et al. [36, 40]. In summary, in that work, it was recognized that the limiting factor for cultivated area estimation in small plot agriculture in Africa is the spatial resolution. A possible way to overcome this limitation is on the synergetic use of sensors with different spatial resolutions and characteristics, therefore by optimizing the spatial and temporal resolution in a way that both dynamics are taken into account (Cf. Estimation of cultivated areas using multi-temporal SAR data).

Table 2.

Confusion matrices forest area product – (top) Multi-year ALOS PALSAR-1 FBD; (bottom) Multi-year ALOS PALSAR-1 FBD and seasonal ENVISAT ASAR.

Acknowledgments

We are grateful to the Japan Aerospace Exploration Agency, the European Space Agency, NASA, and the Italian Space Agency (ASI) for the provision of the ALOS PALSAR-1, ENVISAT ASAR, and Cosmo-SkyMed data. Parts of the work included in this article have been carried out within the JAXA Kyoto and Carbon initiative and the ESA Global Monitoring for Forest Security project. Support is acknowledged from the Gordon and Betty Moore Foundation, The David and Lucile Packard Foundation, and the Google.org Foundation.