Multi-sensor forest vegetation height mapping methods for Tanzania

ABSTRACT This paper proposes a new method for mapping of forest cover in Tanzania, in the form of yearly estimates of average vegetation height from time-series of Landsat and ALOS PALSAR satellite images. By using airborne laser scanning data and Landsat-8 data from 2014, a regression between average vegetation height and the specific leaf area vegetation index is established. By using all available Landsat acquisitions of the same area within 1 year, and producing a yearly estimate of vegetation height, the estimation error variance is reduced. The variance is further reduced by Kalman filtering the sequence of yearly estimates. A multi-sensor version of the method comprises application of the radar backscatter when L-band SAR data is available. To conclude, we have demonstrated that estimation of mean vegetation height is possible from dense time series of optical and SAR satellite data. Change detection was able to detect areas with total loss of biomass.


Introduction
The current rate of deforestation in tropical regions needs to be reduced in order to decrease CO 2 emissions and preserve biodiversity. Advances in remote sensing methods are needed in order to accurately estimate the state of the forests and how they change over time. For the purpose of tracking biodiversity, ten variables have been suggested by Skidmore et al. (2015): (1) species occurrence, (2) plant traits (including specific leaf area and leaf nitrogen content), (3) ecosystem distribution, (4) fragmentation and heterogeneity, (5) land cover, (6) vegetation height, (7) fire occurrence, (8) vegetation phenology (variability), (9) primary productivity and leaf area index, and (10) inundation.
The Landsat archive provides optical multispectral images in 30 m pixel spacing from 1985 and into the near future. Sentinel-2 was launched in 2015 and provides optical multispectral images in 10 and 20 m pixel spacing and with the same wavelength bands as Landsat-8. Therefore, methods that can estimate essential forest parameters from Landsat and Sentinel-2 are important in an operational system for the measuring, reporting and verification (MRV) of forests in tropical countries, including Tanzania, where the current research was conducted.
It is not immediately clear that vegetation height may be estimated from Landsat images. If a pixel is green, it is difficult to estimate whether the land cover is grass, tall trees or agricultural crops. However, given the seasonal variations in rainy and dry periods that dominate in Tanzania, one may assume that trees are better at preserving their greenness than grass and crops. So, by observing the vegetation repeatedly during 1 year, one may be able to estimate the vegetation height from a dense time series of multispectral optical satellite images such as Landsat and Sentinel-2. Hansen and Loveland (2012) points out that the recent opening of free access to the Landsat archive, coupled with advances in computing power, means that change analysis will no longer be based on comparison of individual scenes. Rather, automated image pre-processing and land cover characterization methods are needed. Hansen et al. (2013) generate global annual maps of per cent tree cover and forest gain and loss at 30 m pixel spacing from Landsat data for 2000-2012. Zhu, Woodcock, and Olofsson (2012) use all available Landsat images to monitor forest disturbance over a 2-year period in a study area in South Carolina and Georgia, USA. The seasonal variation during 1 year is modelled and used to predict the pixel state in the current image, presuming no change in land cover. The predicted value is compared with the actual value to determine if a change had occurred. Three repeated observations are required to confirm a change. However, this method may fail CONTACT Øivind Due Trier trier@nr.no Section for Earth Observation, Norwegian Computing Center, Gaustadalléen 23A, P.O. Box 114, Oslo NO-0314, Norway in Tanzania for two reasons. First, in order to model the seasonal variations, the method requires a land cover classification, but this was not available to us. Second, the observed variations in Tanzania are sometimes deviating from a clear seasonal pattern (Figure 1). By comparing four Landsat images of the same area in Tanzania, all acquired in the same month but in different years, one may falsely assume that large areas have been deforested from 1987 (Figure 1(a)) to 1997 (Figure 1(b)). However, only 1 year later (Figure 1(d)), the vegetation suddenly looks more vigorous than in 1987. Then, after another 5 years (Figure 1(d)), the vegetation looks less vigorous than in 1987 and 1998, but more vigorous than in 1997. The reason may be due to varying amounts of rain, or, in occasional years, even lack of rain, in the so-called short rain period, which usually takes place from late November until January.
In the present study, the focus was on monitoring the tree vegetation. Although optical images are useful for monitoring the presence of, e.g. plant chlorophyll by computing for example the normalized difference vegetation index (NDVI; see e.g. Aase & Siddoway, 1981;Campbell, 2006;Townshend, Goff, & Tucker, 1985), there is no reliable way to extract forest variables like the forest vegetation height or the biomass directly from a single Landsat image of an area in Tanzania. One strategy could be to use several repeated acquisitions during the same year, and look for seasonal change patterns in, e.g. greenness, dryness, or some vegetation index, in order to learn a statistical relation between seasonal changes and a forest variable of interest. As hundreds of repeated 30 m pixel spacing Landsat acquisitions from 1984 to date are available for all of Tanzania, it seems obvious that this data source must be exploited, possibly in combination with other data sources.
Synthetic aperture radar (SAR) images overcome cloud cover problems, since the radar wavelengths penetrate almost all clouds. SAR images of Tanzania have been acquired in three different wavelengths: X-band, C-band and L-band. Of these, L-band SAR penetrates deepest into the forest canopy and is therefore best suited to estimate biomass (Mitchard et al., 2009) by capturing large tree structures like stems and thick branches which constitute the dominant fractions of the total tree biomass, ignoring low herbaceous vegetation, thin branches and foliage.
Estimation of vegetation height may be done with a variety of remote sensing technologies (Table 1). Airborne laser scanning (ALS) is currently the most accurate method. The laser scanner records the travel  time from a pulse is emitted until one or more reflections return to the sensor. With accurate recording of the aircraft's position and each laser beam's direction during the scanning, the recorded times may be converted to (x, y, z) positions. If a laser pulse hits a tree canopy that is not too dense, then one part of the pulse may be reflected back towards the sensor, while another part of the laser pulse may continue deeper into the canopy and may even hit the terrain surface before being reflected. A number of reflections from the canopy interior may also occur. Naesset (1997) uses ALS to estimate mean tree height at the stand level, with coefficient of determination R 2 = 0.94 and root mean square error (RMSE) = 7.7%. The sampling is very sparse, with 15 cm footprint diameter and 3 m average distance between footprints, thus only 0.2% of the area is sampled. However, with improved laser scanner technology, sub-metre sampling distance has been available for more than 10 years. Kwak, Lee, Lee, Biging, and Gong (2007) use ALS to detect individual trees and estimate their heights, with R 2 = 0.77 and RMSE = 8.2%. Hopkinson, Chasmer, Lim, Treitz, and Creed (2006) use an indirect method to estimate canopy height from ALS data: a regression on the standard deviation of the difference between first and last returns. The R 2 performance is very good (0.95) but the root mean square error exceeds 10% of the measured mean height.
Many authors (e.g., Ioki et al., 2014;Kwak et al., 2007;Lefsky, Keller, Pang, Camargo, & Hunter, 2007) use ALS data to estimate a canopy height model. First, the individual ALS points are classified as one of "ground", "vegetation", and, optionally, other classes like "building", etc. This classification may be done using the TerraScan software (http://www. terrasolid.com). The ALS points labelled as "ground" are used to estimate a digital terrain model (DTM). The ALS first return points of both "vegetation" and "ground" are used to estimate a digital surface model (DSM). A canopy height model is obtained by subtracting the DTM from the DSM.
Although ALS is a direct measurement of terrain and vegetation height, the resulting canopy height model might underestimate the true canopy height (Lefsky et al., 2007). First, since ALS is a sampling technique, the highest locations on individual tree canopies may sometimes be missed, especially with cone-shaped canopies of, e.g. spruce. Second, the terrain surface may be covered by dense understorey vegetation, thus precluding the terrain surface from being sampled, so that the understorey vegetation may be mistaken as the terrain surface.
The main limitation of ALS is that it is cost prohibitive to be applied wall-to-wall for larger territories like the entire territory of Tanzania (945 000 km 2 ), not to mention the need of repeated acquisitions to monitor change. Thus, alternative remote sensing methods need to be investigated. Véga and St-Onge (2008) use stereo-pairs of aerial images to map historical canopy height models. Image matching techniques are used to estimate DSMs for each of the years 1945, 1965, 1983 and 2003. A recent ALS survey is used to estimate a common DTM. The accuracy is close to that of ALS, and has a huge potential in creating historical records of forest canopy height, given the vast amount of stereo aerial photography in many countries. Lefsky et al. (2007) use the Geoscience Laser Altimeter System (GLAS) sensor onboard the Ice, Cloud and Land Elevation Satellite (ICESat) to estimate forest canopy height. The method does not need an external DTM. The RMSE = 25% is not as good as with ALS. Another limitation of the method is that the sampling is only along transects. To overcome this, Lefsky (2010) combined ICESat GLAS with Modis 500 m data to produce a global forest canopy height map, but with higher RMSE = 43%. Zhang et al. (2014) also use ICESat GLAS, but in combination with Landsat. They argue that there is a linear relationship between leaf area index (LAI) and forest canopy height. Thus, by estimating LAI from Landsat data and calibrating with GLAS point observations of canopy height, a continuous forest canopy height estimate is obtained at 30 m resolution, with RMSE = 35%. Kellndorfer et al. (2004) extracts a DSM from the shuttle radar topography mission (SRTM), which made a global DSM at 90 m grid cell resolution using C-band SAR in interferometric mode. By using the national elevation data of the United States as the DTM, a canopy height model is estimated. Kugler, Schulze, Hajnsek, Pretzsch, and Papathanassiou (2014) investigate the use of TanDEM-X interferometric SAR for forest height estimation in three different forest types. Two identical satellites with X-band SAR instruments share the same orbit. One is emitting a pulse and both are receiving. From this, a DSM is obtained at 12 m resolution. In single polarization mode, an external DTM is needed. In dual polarization mode (HH+VV), an external DTM is not needed, provided some part of the radar pulse is able to penetrate the forest cover, reach the terrain surface and reflect back to the satellites. In both polarization modes, for about 20% of the pixels, the method is not able to provide an estimate. Lei and Siqueira (2014) estimate forest height from ALOS PALSAR L-band SAR, using interferometric correlation magnitude, with R 2 = 0.76 and RMSE = 18.8% for mixed forest in Maine, USA. Sexton, Bax, Siqueira, Swenson, and Hensley (2009) compare using ALS, SRTM, and airborne Xand P-band SAR interferometry (GeoSAR), for forest canopy height estimation in pine and hardwood forests in North Carolina, USA. The ALS data has 0.5-1.0 m diameter footprints spaced at 4-6 m, which is much coarser than in the ALS studies of Kwak et al. (2007) and Hopkinson et al. (2006). In pine forest, ALS performs much better (R 2 = 0.83, RMSE = 7.9%) than GeoSAR (R 2 = 0.70, RMSE = 28%), which again performs much better than SRTM (R 2 = 0.54, RMSE = 48%). In hardwood forest, the performance was poorer than in pine forest, but GeoSAR (R 2 = 0.38, RMSE = 33%) is still better than SRTM (R 2 = 0.35, RMSE = 66%). Helmer et al. (2010) estimate current vegetation height from a time series of images from Landsat (1984Landsat ( -2003 and EO-1 Advanced Land Imager (ALI) (2005)(2006), using a regression model. The independent variables in the regression are the NDMI vegetation index and forest disturbance indicators (undisturbed, goat grazing, cleared for agriculture, and fire disturbance) The R 2 = 0.87 is good, but RMSE = 25% is high. Wilkes et al. (2015) combine spring and summer Landsat images with a 10 year Modis NDVI time series to estimate forest canopy height in forests in eastern Victoria, Australia, with R 2 = 0.58 and RMSE = 31.3%. Better results are obtained by Wolter, Townsend, and Sturtevant (2009), who estimate tree height from two SPOT-5 multispectral images with 10 m resolution. They use regression on vegetation indexes and local statistics, obtaining R 2 = 0.92 and RMSE = 9.4% on a coniferous forest and R 2 = 0.69 and RMSE = 7.9% on a hardwood forest.
Chopping, Nolin, Moisen, Martonchik, and Bull (2009) estimate forest canopy height from the Multiangle Imaging Spectroradiometer (MISR), which uses four forward and four aft viewing angles in addition to straight down. The results (R 2 = 0.66, RMSE = 51%) are the least convincing of the studies in Table 1. Kalacska et al. (2007) estimate canopy height from EO-1 Hyperion satellite data with 220 hyperspectral bands at 30 m grid spacing, acquired over a deciduous forest in the dry season. Regression is done on predictors extracted from wavelet decomposition, and has a better coefficient of determination (R 2 = 0.92) than regression on a single vegetation index (R 2 = 0.66).
Cartus, Kellndorfer, Rombach, and Walker (2012) compare using ALS, PALSAR, Landsat-7, and a combination of PALSAR and Landsat for canopy height estimation in pine plantations in Chile. Using PALSAR data from several dates in a single year (2007), with single (HH) and dual (HH+HV) polarization and interferometric coherence, R 2 = 0.72 and RMSE = 11.5% are obtained. By including three cloud-free Landsat acquisitions from 2005-2007, this is improved to R 2 = 0.83 and RMSE = 10.2%. Using Landsat alone results in R 2 = 0.64 and RMSE = 13.5%. Hyyppä et al. (2000) compare stand-level mean height estimation using different airborne and satellite sensors. Of these, airborne C-band full polarimetric SAR performed best, with R 2 = 0.77 and RMSE = 18%. For the optical data, airborne hyperspectral data with 30 channels and 2 m resolution performed best, R 2 = 0.47 and RMSE = 31%. Despite different resolutions, Spot multispectral satellite images with 20 m resolution and three channels, green, red, and near-infrared, performed equally well (R 2 = 0.37, RMSE = 35%) as scanned aerial photographs in 0.85 m resolution and the same three channels (R 2 = 0.34, RMSE = 33%). However, Landsat 30 m resolution performed slightly worse (R 2 = 0.26, RMSE = 37%) despite having more spectral bands. The ERS-1/2 tandem mission provided interferometric C-band SAR acquisitions with 1 day time difference. From these data, coherence images at eight different date pairs, each with 1 day lag, were used, but gave poor results, R 2 = 0.15 and RMSE = 39%.
To summarize, for ALS, the coefficient of determination, R 2 varies between 0.77 and 0.95 with RMSE between 7.1% and 10.6%; Landsat, R 2 0.26-0.87 with RMSE 13.5-36.7%; Palsar, R 2 0.72-0.76, RMSE 11.5-18.8%; other satellite sensors/combinations, R 2 0.15-0.92, RMSE 6.5-51% (Table 1). The large variations in R 2 and RMSE are in part due to the different sizes of the study sites, which range from a single forest plantation to, e.g. the entire state of California, USA. In many of the studies, a canopy height model from ALS data is used as the true values instead of field data.
In this study, we hypothesized that it should be possible to produce yearly estimates of average vegetation height at 30 m pixel spacing from a time series of Landsat and ALOS PALSAR data. The objective was therefore to develop such a method that could be applied operationally over large areas of tropical forests. For calibration of the method, we use ALS data with 1.8 first returns per square metre on average to estimate a canopy height model at 1 m pixel spacing. When aggregated to 30 m pixel spacing, this is regarded as exact measurements of mean vegetation height.
To reduce noise and to provide estimates for areas that have not been observed due to cloud cover, some sort of smoothing may be needed. Since we do not want coarser resolution than 30 m, smoothing in the spatial domain should be avoided. Two popular methods for smoothing in the time domain is Kalman filtering (Kalman, 1960) and hidden Markov model. Trier (2011, 2012)) have used a hidden Markov model for temporal analysis of tropical rain forest in Brazil. However, by using the same approach in Tanzania, the seasonal fluctuations of the optical signal may falsely be interpreted as state changes in a hidden Markov model. Kleynhans et al. (2011) used extended Kalman filtering on an NDVI time series for land cover change detection.

Study areas
Two study areas were included: Liwale and Amani ( Figure 2). Liwale is representative for 90% of Tanzania, while Amani was included as an extreme case in order to explore the limitations of the proposed method.
The Liwale study area (S9°54ʹ, E37°38ʹ) covers 15,867 km 2 of Miombo woodlands with altitudes in the range 150-900 m above sea level and is located in the Lindi region, southern Tanzania. Miombo woodlands cover 90% of Tanzania. As described by Mauya et al. (2015), the rainfall pattern in Liwale is bi-modal with a dry season from June to October. A short rainy period usually starts in late November and lasts until January. There is a dry spell in February followed by a longer wet season which lasts from March until May. Mean annual rainfall is 800 mm.
The Amani study area (S5°08ʹ, E38°37ʹ) covers 127 km 2 with altitudes in the range 200-1150 m above sea level. It includes the 83.6 km 2 Amani Nature Reserve, which is a tropical submontane rainforest. The Amani area is located in the East Usambara Mountains, which are parts of the Eastern Arc Mountains in eastern Tanzania. In this forest ecosystem, rain falls throughout the year with two wet seasons, April to May and October to November, and the area receives around 2000 mm rainfall per year (Hansen, Gobakken, Bollandsås, Zahabu, & Naesset, 2015).

ALS data
ALS data were collected by TerraTec AS, Norway, in 2012 for both Liwale and Amani, and repeated in 2014 for Liwale (Table 2). The ALS data in Liwale ( Figure 3) were collected along 1.4 km wide and parallel strips in a systematic design at 5 km intervals, i.e. there was a gap of 3.6 km between neighbouring strips. A 113 km (east-west) × 156 km (north-south) area was covered with 34 strips in the east-west direction. The repeat acquisition in 2014 was done along the same strips to allow for change analysis. The ALS data in Amani ( Figure 4) covered 127 km 2 wall-to-wall, including a protected high-biomass montane tropical rain forest with some trees taller than 50 m. This area is not typical for Tanzania, but represents an extreme situation with rain forest. The data were collected on six dates in January-February 2012.
From each ALS data set, average vegetation height was computed on the 30 m Landsat grid. In the ALS data, each ALS pulse had been labelled with class ("ground" or "other") and return number. From all "ground" returns, a  DTM was created at 1 m pixel spacing. From all the first returns, a DSM was created at the same pixel spacing. By subtracting the DTM from the DSM, a normalized DSM (nDSM) was obtained, which may be used as an estimate of vegetation height. By aggregating this to the 30 m pixel spacing Landsat grid, the resulting average vegetation height map is regarded as exact for the purpose of developing a method to estimate average vegetation height from satellite images with 30 m pixel spacing. However, there is a potential problem with the 2014 ALS data in Liwale, since it was acquired in the dry period, which may allow the ALS pulses to penetrate more into the canopy before the first returns occur. Field work conducted at the times of the ALS campaigns (Ene et al., 2017) gave an estimated biomass loss in the range 0.11-0.35Mg/ha from 2012 to 2014, depending on which estimator was used. None of these changes are statistically significantly different from 0 at the 95% confidence level. Therefore, we will scale the 2014 nDSM by a factor (Table 3)  Only Landsat data that were terrain corrected (level 1 terrain corrected, L1T) were used; thus, a number of acquisitions that were geocoded but not terrain corrected (level 1 geocoded, L1G) were discarded from the analysis (Table 4)  2007-2010 with fine beam dual polarization data (FBD) in HH and HV polarization with satellite path numbers 566, 567 and 568. The level-1.1 data have been processed, i.e. geocoded, radiometrically corrected, including slope correction, with Norut's in-house developed SAR pre-processing system GSAR (Larsen, Engen, Lauknes, Malnes, & Høgda, 2005) into the same 30 m pixel size grid as the Landsat data.

Methods
The overall goal of the proposed method is to estimate forest vegetation height all over Tanzania, despite obvious challenges due to the following: (1) Varying cloud cover. It is very rare that a single Landsat image is entirely cloud free. (2) Seasonal variations. The yearly cycle usually has two rain periods and two dry periods. However, these vary in length, intensity and spatial distribution from year to year.
(3) Saturation. The method may underestimate the vegetation height when the true average vegetation height is above some limit.
The approach to arrive at a general method is as follows: (1) Establish a method that works for the Landsat data that overlaps the ALS area in Liwale, for the same year (2014).
(2) Modify the method so that it may be used for Landsat data that overlaps the ALS area in Liwale, from any year. (3) Extend the method so that it may be used for Landsat data everywhere in Tanzania.
The Amani area is not used for method calibration, but is included in the evaluation as an odd case.  The method is novel in the following aspects.
(1) It uses all available Landsat acquisitions. (2) It provides vegetation height estimates for all years in the time series.

Vegetation height from Landsat images
There is a linear relationship between LAI and forest canopy height, provided the forest canopy is not closed or nearly closed (Zhang et al., 2014). We used the specific leaf area vegetation index, SLAVI, (Lymburner et al., 2000) as an estimate of LAI. SLAVI uses three Landsat bands. The index uses the positive correlation of leaf chlorophyll with the NIR band, and the negative correlation with the red band. In addition, the short wave radiation is absorbed by leaf water content, and the effect is stronger in the SWIR 2 band (2.2 µm) than in the SWIR 1 band (1.6 µm). SLAVI is calculated by Compared with other vegetation indices like NDVI and NDMI (Jin & Sader, 2005), SLAVI is less affected by cloud shadows, which may remain even after cloud and cloud shadow masking.
In order to establish a relationship between SLAVI and forest canopy height, the nDSM with 30 m pixel spacing from the 2014 ALS data from Liwale were used together with SLAVI maps from all 21 Landsat-8 acquisitions from 2014 of the same area. The normalized DSM (nDSM) provides an estimate of the vegetation height. The nDSM was aggregated to the 30 m Landsat grid, thus providing an average vegetation height within each Landsat pixel. After the scaling of the 2014 nDSM, as described above, it was regarded as an exact measurement of mean vegetation height for each Landsat pixel.
In order to investigate if there is a relationship, which may be used in regression, between SLAVI and vegetation height, the following plot was made. The SLAVI axis was divided into bins of width 0.1. Then, for each Landsat-8 image from 2014, only the cloud-free pixels that had a mean vegetation estimate from the 2014 nDSM of Liwale were used. Of these, for all pixels within one SLAVI bin, the average of the nDSM vegetation height estimates for those pixels was used. This was repeated for all SLAVI bins. This resulted in the 21 curves in Figure 6, with one curve for each Landsat acquisition date.
For each curve (Figure 6), there appeared to be a near-linear relation up to a peak. Beyond the peak, the relation seemed to be near-linear with a negative slope. However, by closer inspection of the scatter plots of the data (Figure 7), there was only a small fraction of data points beyond the peak. Thus, this part of the curve may be ignored. Scatterplots of SLAVI versus vegetation height for each of the 21 Landsat images (e.g. Figure 7) also indicated that there may be a linear relation between SLAVI and vegetation height for most of the images, but also that the variance is high. For each Landsat acquisition during 2014, a regression between SLAVI (based on Landsat-8 bands 4, 5 and 7) and the acquired ALS vegetation height was established (Table 5, columns 1-5): where y is vegetation height, x is SLAVI, a is the y axis intercept at x = 0 (a theoretical vegetation height for SLAVI = 0) and b is the slope of the linear relationship. However, when we encounter Landsat images for other years, with no ALS data available, it is difficult to estimate both a and b. It would be better to have a fixed a = a 0 , so that only one parameter, b, needs to be estimated for each new image. This was achieved by the simultaneous regression  (Table 5, columns 8-9) gives a common regression constant a 0 = −3.40. There is no direct physical interpretation of this value; a negative vegetation height is not possible in reality. For each of the N = 21 Landsat-8 images of 2014, the two regression lines may be compared visually with the scatter plot ( Figure 8). The simultaneous regression gives a slight increase in root mean square error (Table 5, column 7 vs. 11).
For a new image k, from another year than 2014, we do not perform a regression between an ALS canopy height model and the SLAVI values, as the true vegetation height values may deviate substantially from the 2014 ALS estimate, especially if there is a substantial number of years between the image acquisition and the ALS acquisition. Instead, we wanted to use the regression Here, the regression slope b k is unknown. However, by comparing the new image k to a reference image r, the regression slope was estimated as where b r is the regression slope of a reference image of the same Landsat path/row, and µ r and µ k are the mean values of the simultaneously cloud-free pixels of the reference image and the new image, respectively. For example, if the new image k has a mean SLAVI value of only half of that of the reference image r, then the regression slope b k of the new image k will be twice the slope b r of the reference image. The effect is a data scaling, but with a common offset (regression constant). An alternative approach could have been to use a common regression slope and change the regression constant for each new image. This would have been a data shift, i.e. a k = a r + (µ r -µ k ).  The reference image must be almost cloud free to allow for a large number of simultaneously cloud free pixels. To establish a new reference image for another Landsat path/row, the corresponding regression slope was estimated using the same estimator, but with the mean values estimated for all cloud free pixels in each image. By cloud free, we mean that Fmask has labelled them as valid, and not any of: cloud, cloud shadow, water region, or outside image.
For each image of the same Landsat path/row within 1 year, the regression slope and then the vegetation height for all cloud free pixels were estimated. The yearly average for each pixel was obtained by computing the mean value of the height estimates for all cloud free observations of that pixel within the year. The number of cloud free observations of each pixel was stored in a separate image for later use by the Kalman filtering.
The Kalman filtering requires estimates of the variance associated with the yearly average vegetation height estimates. For 2014, we have measurements of average vegetation height, µ i , for pixels i = 1, . . ., I. For each pixel i, we have N i cloud-free satellite measurements y ij in the same year, j = 1, . . . N i . The yearly average of the average vegetation height for pixel i is The N i measurements y ij of pixel i are not independent. We model this as where µ i is the true value, ε i is the dependent part of the measurement error and ω ij is the independent part of the measurement error. We will use σ 2 and τ 2 to denote the error variances: For the yearly estimate y i for pixel i, we then have This means that we may reduce the variance of y i by adding more observations (i.e. increasing N i ), but also that we will never get rid of σ 2 . However, this estimate of σ 2 is a worst-case estimate. To estimate σ 2 and τ 2 we usê N i is the total number of observations and I is the number of pixels for which we have true values μ i . For the 2014 data, this givesσ 2 þτ 2 ¼ 10:26 and a worst-case estimate ofσ 2 ¼ 8:05, which giveŝ τ 2 ¼ 2:21 (Table 6). This means that a large number of observations within 1 year has a limited effect on reducing the estimated variance given byσ 2 þ 1 N iτ 2 . However, for the purpose of Kalman filtering, the estimated variance (per pixel) is used to indicate to what extent the new observations are to be trusted. We want to differentiate more between few and many observations, and do this by reducing the estimate for σ 2 and increasing the estimate for τ 2 . We suggest to usê σ 2 ¼τ 2 ¼ 5:13, which may better reflect the added value of having more observations within 1 year.

Vegetation height from PALSAR images
Vegetation height was also estimated for L-band PALSAR data for the years 2007-2010. SAR volume scattering increases with vegetation and it has been shown that there is an exponential-like relationship between SAR backscatter and aboveground biomass (Mitchard et al., 2009). It is therefore reasonable to assume that there is also an exponential-like relationship between backscatter and vegetation height. In order to calibrate this relationship in our data, the ALS data were aggregated to the pixels of the ALOS PALSAR mosaic. In order to reduce for speckle effect, the average ALOS PALSAR mosaic over the 2007-2010 period was aggregated into 90m resolution pixels.
For HH polarization the regression model was h vegt;HH ¼ 3557:7e 0:4197σ HH with R 2 = 0.9408. For HV polarization, the regression model was h vegt;HV ¼ 944:09e 0:5283σ HV with R 2 = 0.9420. These regression models were used to estimate the vegetation height for each year 2007-2010 individually in 90 m pixel spacing. The results based on HH and HV polarization individually were then averaged for a final combined result. Finally, it was resampled in 30 m pixel We observe that there was a declining vegetation height trend through the 4 years. As the 2010 vegetation height mean value is reasonably close to the ALSderived mean value, and closest in time to the 2012 ALS data, it seems reasonable to use the 2010 PALSAR data in order to estimate the varianceσ 2 P of the mean vegetation height predictions from PALSAR. σ 2 P ¼ 7:63m 2 This variance will be used as input to the Kalman filtering (below). The Kalman filtering also requires a gain factor. As the estimated Kalman gain factor for 2010 was close to 1.0 (Table 7), we used the value 1.0 in the subsequent analysis. This was also in agreement with the declining trend of the average vegetation height, i.e. the mean of the vegetation height was expected to be slightly reduced from 2010 to 2012.

Reducing noise
The vegetation height estimates from individual years may vary substantially, given the varying number of satellite observations within each year, and the varying cloud cover in Landsat images. In order to reduce fluctuations in the estimated vegetation height that are not caused by actual changes in vegetation height due to regrowth/growth and anthropogenic or natural losses of trees and other vegetation, but related to biophysical response of the vegetation to the variability in climatic conditions, Kalman filtering was applied. We assume no change as the normal situation, i.e. regrowth/growth is balanced by anthropogenic or natural losses of trees and other vegetation.
The Kalman filtering was designed so that it may accept any of the following four situations for each year in the sequence: • No data • Estimate from Landsat • Estimate from PALSAR • Estimate from both Landsat and PALSAR The Kalman filter was used to model a sequence of process states. The process state y k+1 at time t k+1 is dependent on the state y k at a previous time t k : Figure 9. PALSAR backscatter versus ALS vegetation height. Exponential regression lines (solid lines) were fitted to the data points (cyan symbols for HV and purple symbols for HH polarization). where F is the change factor and w k~N (0,Q) is the process noise, assumed to be normally distributed with zero mean value and variance Q. The observation x k at time t k is: Here, H is the gain factor for the measurement, and v k is the measurement error. We have two measurements, x P,k from PALSAR and x L,k from Landsat: So, for each pixel, we measured the yearly values of x L,k from Landsat and/or x P,k from PALSAR for all the years k for which we have satellite observations, and estimated the sequence of yearly y k vegetation height. The following parameters were used: Of these, R L and R H were estimated from the data, as described above. F = 1.0 was used to indicate that no trend is assumed, i.e. if a pixel has no cloud-free observations, then the estimated mean vegetation height will remain unchanged. Still, using F = 1.0 does not preclude changes. Q may be estimated by comparing the 2012 and 2014 nDSMs. The 2-year change may be modelled using Since F = 1.0, Since w k has 0 mean value, the variance of w k þ w kþ1 is 2 times the variance Q of w k : This gave 2 × Q = 1.548, i.e., Q = 0.774. However, in the Kalman filtering, Q = 0.774 cannot be used, since such a low value does not allow for total loss of biomass (e.g. deforestation) to happen. The low value may be optimal on average, since total loss only happens to a small fraction of the pixels, yet total loss in the form of deforestation is an important event to map. By trial and error, we found that Q = 10 was needed to capture events of total loss of biomass.

Change mapping
The output of the Kalman filtering was a smoothed version of the time series, with yearly estimates of vegetation height. Change maps between any 2 years in the time series may be produced by computing the difference. However, one should avoid including years near the beginning of the time series, since only a limited amount of smoothing has been done.

Results
The deviation of the mean is low for all three versions of the method (Table 8). For change mapping, the method was able to capture most of the events of total biomass loss ( Figure 11). For the Amani area (Figure 4), which is a montane tropical forest with some trees taller than 50 m, the proposed method performs less well. PALSAR data could not be used for vegetation height estimation due to the steep topography. Therefore, only Landsat data were used: path 166, row 064 from 1984-2012. On average, for Landsat 2012 data, the method underestimated the mean vegetation height by 6.6 m ( Table 9). Kalman filtering did not improve the bias, and the variance also remained high. However, Kalman filtering was able to provide estimates for 6814 pixels that lacked 2012 Landsat-7 coverage (Figure 12(a,b); Table 9). From the difference image (Figure 12(d)), we observe that the method overestimated the vegetation height in areas with low vegetation. For areas with high vegetation, the method underestimated the average vegetation height by up to 40 m. This is likely a saturation effect.

Discussion
The results of average vegetation height estimation show that it is possible to extract height information from Landsat images of forested areas provided the forest is not too dense. By using time series instead of individual images, the variance is reduced substantially. By adding PALSAR data to the time series, the variance might be further reduced.
For change detection, the results are less convincing. Areas with total loss of biomass were detected, but gradual changes are not captured well. The difficulties are due to the variation of the yearly height estimates from satellite. This variation is often larger than the actual change in vegetation height as measured by the ALS data.
The performance of the proposed method in terms of root means square error (Tables 8 and 9) is not convincing compared to existing methods (Table 1). However, most of these studies are on forests only, while our method has been applied on all land cover types, including agricultural fields and mixes of agriculture and tree vegetation.   We have based our study on a simple, linear regression between one vegetation index (SLAVI) derived from three satellite bands, and mean vegetation height. A regression on multiple vegetation indexes could be used with the goal to reduce the variance.
The regression between PALSAR mean backscatter and ALS (Figure 9) seem to deviate systematically from the fitted exponential curve, in that high vegetation values >17 m are underestimated while medium vegetation values 3-12 m are underestimated. Thus, a different parametric nonlinear function may produce better results.
On the other hand, a potential problem with the regression between PALSAR and ALS is that the PALSAR data are from 2007 to 2010 while the ALS data are from 2012. Thus, vegetation height changes between the PALSAR and ALS acquisitions are ignored. These changes may be a mix of height increase due to growth and decrease due to felling of individual trees and clear cuts. If we assume that high trees are more likely to be felled than small trees, while average height growth is more likely to occur in areas with low average vegetation height, then the systematic deviations observed in Figure 9 may be explained by the time difference of acquisitions.
Another way of reducing the variance is through more averaging, in time and/or space. Given the low number of acquisitions per year in 1984-1990, one may use the estimate for 1990 after Kalman filtering as the first year in the time series. For the purpose of reporting of carbon stocks, Tanzania needs to provide one single number every second year for the entire country. Since 90% of Tanzania is Miombo woodlands, on may use a regression between biomass and average vegetation height developed for Miombo woodlands. For each reporting year, one may run the processing chain on all the 48 path/row pairs that are needed to Cover Tanzania, mask areas outside of Tanzania, and produce one estimate for each path/row and then aggregate again to produce one single estimate for the entire country.
One strategy to obtain more optical observations from the past than may be provided by Landsat could be to include daily acquisitions at a coarser resolution, e.g. Modis (250 m-1 km resolution), VIIRS (375-750 m) or Sentinel-3 SLSTR (500-1000 m) and OLCI (300 m). It is straightforward to modify the Kalman filtering step to include more sensors. However, all images must be resampled to the same resolution.
To improve the accuracy of the method for present and future mean vegetation estimates, it seems obvious to include Sentinel-2 data, with 10 and 20 m pixel spacing instead of 30 m, and more frequent revisit. Another problem that Sentinel-2 may help to solve is that two 30 m × 30 m Landsat pixels with the same average height may have different spatial distributions of vegetation height. An area with 50% canopy cover of 20 m tall trees may have a different spectral and temporal signature than an area with 100% canopy cover and 10 m tall trees. By dividing each Landsat pixel into 3 × 3 Sentinel-2 pixels, these differences may become apparent, unless the 20 m high tree canopies that occupy 50% of one Landsat pixel are homogeneously distributed over the 3 × 3 Sentinel-2 pixels.
The method is developed for open forests and woodland landscapes, with gradual changes in tree crown cover. The method was calibrated with ALS data from Liwale, since this landscape is representative for the majority of forests and woodland landscapes in Tanzania. The method saturates at around 10 m average vegetation height. For tropical rain forests, like in Amani, large areas had too low estimates. However, the method may still work for the detection of logging, of which many may be illegal conducts when it takes place in rain forests in Tanzania.
Accurate estimation of vegetation height and biomass for all of Tanzania by ALS is not considered as a realistic alternative, given the high acquisition cost. Forest change detection is also possible with TerraSAR Tandem-X, but again the cost is prohibitive at the moment. Given the free access to optical satellite data at various resolutions, further research into automatic methods for the extraction of essential forest variables seems to be a viable alternative.

Conclusions
To conclude, we have demonstrated that estimation of mean vegetation height is possible from dense time series of optical and SAR satellite data. However, better accuracy in the height estimates is desired. With better resolution and higher acquisition frequency, Sentinel-2 may provide exactly that. We recommend further studies on applying the developed method on Sentinel-2 data.