Mapping Forest Canopy Height over Continental China Using Multi-Source Remote Sensing Data

Spatially-detailed forest height data are useful to monitor local, regional and global carbon cycle. LiDAR remote sensing can measure three-dimensional forest features but generating spatially-contiguous forest height maps at a large scale (e.g., continental and global) is problematic because existing LiDAR instruments are still data-limited and expensive. This paper proposes a new approach based on an artificial neural network (ANN) for modeling of forest canopy heights over the China continent. Our model ingests spaceborne LiDAR metrics and multiple geospatial predictors including climatic variables (temperature and precipitation), forest type, tree cover percent and land surface reflectance. OPEN ACCESS Remote Sens. 2015, 7 8437 The spaceborne LiDAR instrument used in the study is the Geoscience Laser Altimeter System (GLAS), which can provide within-footprint forest canopy heights. The ANN was trained with pairs between spatially discrete LiDAR metrics and full gridded geo-predictors. This generates valid conjugations to predict heights over the China continent. The ANN modeled heights were evaluated with three different reference data. First, field measured tree heights from three experiment sites were used to validate the ANN model predictions. The observed tree heights at the site-scale agreed well with the modeled forest heights (R = 0.827, and RMSE = 4.15 m). Second, spatially discrete GLAS observations and a continuous map from the interpolation of GLAS-derived tree heights were separately used to evaluate the ANN model. We obtained R of 0.725 and RMSE of 7.86 m and R of 0.759 and RMSE of 8.85 m, respectively. Further, inter-comparisons were also performed with two existing forest height maps. Our model granted a moderate agreement with the existing satellite-based forest height maps (R = 0.738, and RMSE = 7.65 m (R2 = 0.52, and RMSE = 8.99 m). Our results showed that the ANN model developed in this paper is capable of estimating forest heights over the China continent with a satisfactory accuracy. Forth coming research on our model will focus on extending the model to the estimation of woody biomass.


Introduction
Forests play a key role in the global climate system and carbon cycle.Forests store carbon in their above-and below-ground biomass [1].As an important predictor of forest biomass and carbon stock, the vertical structure of forests has been well monitored in previous studies [2][3][4].Light detection and ranging (LiDAR) remote sensing, is useful in the large-scale investigations of forest structural attributes, such as forest canopy height [2,[5][6][7][8][9][10].
In recent years, LiDAR instruments have demonstrated their capability to estimate forest heights [11].For instance, the Geoscience Laser Altimeter System (GLAS) on board the Ice, Cloud and Land Elevation Satellite (ICESat) has provided within-footprint height measures at the global scale [2,5,12].However, the lack of valid GLAS data in some regions is problematic.Because the GLAS data are spatially incomplete, multiple geospatial predictors were supplemented to obtain spatially-continuous forest height maps in recently-published approaches [2,3,5,13].In addition, the GLAS data are highly sensitive to topographic features due to its large footprint size, which causes overestimations of forest height [14].In this paper, in order to solve those two problems degrading the accuracy of GLAS-based forest height maps, we proposed a new approach based on the artificial neural network (ANN).
The principal method was to combine GLAS derived tree heights with ancillary variables in the ANN model and then to produce a wall-to-wall forest height map over the China continent.Selected geospatial predictors were a series of climate variables, elevation, vegetation cover and multispectral reflectance data from the Moderate Resolution Imaging Spectroradiometer (MODIS).We resampled all of the input predictors to obtain the final output map of forest heights at a 1-km spatial resolution.

ICESat Data and Processing
The ICESat/GLAS is the first spaceborne LiDAR system, which was designed to obtain characteristics of the Earth's surface structures with unprecedented accuracy [15].The ICESat/GLAS data can provide information related to land surface elevation with a spatial resolution of 70 m (ellipsoidal footprints) and 170 m spaced intervals [16].The global ICESat/GLAS data are available for a period spanning from 2003 to 2008 [17].Due to the short lifetime of GLAS, altimetry information was collected only in the following months: February-March, May-June, and October-November [11,18].There are totally 15 GLAS data products including Level-1 and Level-2 that were disseminated by the NSIDC (National Snow and Ice Data Center).The GLA14 used in this study is an GLAS Level-2 Land Surface Altimetry product that provides information on surface elevations, laser footprint geolocation, waveform parameters and reflectance, as well as geodetic, instrument, and atmospheric corrections for range measurements.In previous studies, the GLA14 product has been used to estimate forest canopy heights within each footprint [1,19].
In this study, we obtained the GLAS data recorded from May to October (2003)(2004)(2005)(2006), because these data represent the approximated growing season and thus, the best leaf-on condition of forests.In our analyses, GLAS footprints over non-forest area were excluded.We further screened invalid GLAS data using several preprocessing filters in order to reduce the impact of slope gradient, atmospheric forward scattering, signal saturation and cloud contamination [14,20,21].The GLAS data were considered valid when footprints were located over the forest classes in the land cover (LC) map with >50% tree cover percent in the vegetation continuous field (VCF).Then, the quality of GLAS shots over the forested area was investigated using GLA14 waveform parameters regarding the cloud, signal saturation and atmospheric forward scattering [1,21,22].Firstly, we removed the invalid GLAS footprints affected by signal saturation and atmospheric forward scattering based on the internal flags ("FRir_qaFlag = 15" and "satNdx = 0") GLA14 product contained [20][21][22].To further remove waveforms with signal dominated by low cloud, we just selected the waveforms based on the filter that was built by using the absolute difference (50 m) between ASTER DEM and the internal elevation value from GLA14 product [21].The terrain slope gradient also impacts the quality of the GLAS data [23].In order to minimize the bias resulting from the terrain slope, we filtered the GLAS shots where the terrain slope is higher than 10 degrees [24].

Land Surface Reflectance
The land surface reflectance data used in this study were acquired from the MCD43B4 which is a product of the MODIS nadir bidirectional reflectance distribution function adjusted reflectance (NBAR) with a spatial resolution of 1 km [25].The NBAR product represents the best characterization of surface reflectance over a 16-day period.In the study, seven spectral bands in the MCD43B4 averaged over the approximated growing season (June-September) between 2004 and 2006 were extracted at a 1-km spatial resolution.

Climate Data
Weather station-based climatic variables, precipitation and temperature, were derived from China Meteorological Data Sharing Service System (CMDSSS).There are total of 754 weather stations that can provide annual meteorological records.The precipitation and temperature data covered a time period from 2004 to 2006.For each weather station, 3-year averaged precipitation and temperature were obtained.Then, we interpolated these two climate variables using ordinary kriging to generate the gridded maps of annual average precipitation and temperature [26].

Ancillary Data
Ancillary data applied in this study involved the digital elevation model (DEM), LC, and VCF.The DEM data were obtained from the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global DEM (GDEM V2) of which spatial resolution is 30 m [27].The use of DEM data is first to remove invalid GLAS shots with a signal contaminated by low clouds and second to generate a gridded map of slope at the same spatial resolution of the GDEM V2 data (30 m).This continuous field of slope was used to correct the GLAS derived tree heights based on the terrain correction method [21].
The LC data were derived from a MODIS global LC product (MCD12Q1, 500-m grid).Additionally, the VCF data were also obtained from the MODIS (MOD44B, 250-m grid).We selected the International Geosphere-Biosphere Programme (IGBP) of the LC and tree cover percent of the VCF over China for the year 2005 [28].

Field-Measured Tree Heights
Field-measured tree heights were used to evaluate our modeled tree heights.We chose 92 plots from three experimental sites.These three field regions are located in Gansu (northwestern China), JiangXi (southern China) and Yunnan (southwestern China), respectively.
The location of the filed experimental regions and corresponding plots is shown in Figure 1.The detailed information on the field sites is delivered in Table 1.Plots in the field sites were selected according to the definition of forest area using LC and VCF in this study.The field measured tree heights were averaged for each plot.

GLAS Tree Height Estimation
The GLAS tree heights were used as training data to build the ANN model for tree height estimation and to validate model performances.GLAS height metrics were extracted based on the "decomposition of GLAS waveforms into multiple Gaussian distribution curves" as described in previous studies [18,19,30].The direct method was used to estimate canopy height based on the vertical difference between the signal start point and the ground peak.The signal start point was provided in the GLA14 product (i_SigBegOff).Also, the GLA14 product contains up to six Gaussian peaks for each GLAS waveform (i_gpCntRngOff).The distance between "i_SigBegOff" and "i_gpCntRngOff" represents the waveform metric RH100 (100% energy return height) of GLAS-derived tree height estimation [9,17].Because the RH100 can be potentially affected by the topographic gradient [20], we applied the topographic correction approach [14,21] to obtain more accurate GLAS tree heights.The calculation formula of GLAS-derived tree heights based on topographic correction is shown as Equation (1), * tan 2 where represents the GLAS derived tree heights, and separately represent the values of the beginning signal and the ground peak of the GLAS full-waveform, d is the GLAS footprint size of 70 m [3,16], and is the topographic slope.

Tree Height Modeling
We employed the ANN algorithm to combine the GLAS-derived tree heights with the geospatial predictors.The neural network algorithms are products of artificial intelligence as black-box models, and the first prototype neural network was proposed in 1943 [31].To date, more than 30 different neural network models have been developed [32], which have been widely used in various fields [33].
Here, we used the feed-forward neural network (FFNN) algorithm for the forest height predictions.The selected model consisted of 11 neurons in the input layer, 11 neurons in the hidden layer and 1 neuron in the output layer.The 11 parameters in the input layer include the LC class, VCF tree cover percent, temperature, precipitation, and seven MODIS NBAR bands.All of the 11 parameters in the input layer are related to forest growth, distribution or the characteristics' expression.The neuron in the output layer refers to the modeled forest canopy heights.Figure 2   To train the ANN model, we applied the back-propagation (BP) process algorithm to train the neural networks [34].For each pixel, we trained the FFNN to estimate forest canopy heights.Based on the temperature and precipitation values, 15 pairs of training data were chosen to train the FFNN while 5 validation pairs remained to prevent the over-fitting of the FFNN.The selection discipline of training and validation data pairs is to obtain 20 GLAS shots that are located in the most similar climatic condition factor (temperature and precipitation) based on Equation ( 2), (2) where is the distance of the climate variable difference between the reference pixel of the canopy height estimation map and the i-th GLAS shot; is the temperature value of the reference pixel; is the temperature value of i-th GLAS shot; is the precipitation value of the reference pixel; is the precipitation value of the i-th GLAS shot; and are respectively the maximum values of temperature and precipitation values over continental China, and and are respectively the minimum values of temperature and precipitation values over the continental China.
After selecting the training and validation data pairs for a pixel, the network gets trained by the training data until the FFNN model loses its best performance given validation data.Then, the trained FFNN model was used for estimating the forest canopy height over the target pixel.

Error Analysis
To assess the performance of the canopy height model in the study, we prepared several reference datasets.These datasets mainly include field-measured tree heights, GLAS-derived tree heights and existing forest height products.The root mean square error (RMSE) was calculated using the following Equation (3): where, is the predicted heights derived from the model; and is the validation heights derived from the reference datasets which were used to validate the predicted tree heights.

Calibration and Comparison with Existing Canopy Height Products
In this study, we present three forms of validation.One is the field calibration of the modeled tree heights, which was based on the field measurements mentioned in Section 2.1.5.For each field plot, we averaged all tree heights to obtain the representative height of this plot.The pixel on the modeled height map that is closest to the field plot centroid was chosen for the validation.The second is that the modeled heights were directly evaluated with the footprint-level GLAS heights (footprint vs. modeled pixel).The third is wall-to-wall map validation from the interpolated GLAS tree height map.Here, we averaged the GLAS tree heights from 20 GLAS shots selected according to Equation (2) to produce a full gridded map.Then, pixel-to-pixel comparisons (FFNN vs. averaged GLAS tree heights) were performed.
Lastly, we compared our modeled tree height map with two existing forest height products.These two reference products are named HSimrad and HNi, respectively, following their creators Simard and Ni [2,3].Our modeled data and two reference products commonly provide forest canopy heights over the continental China, but their modeling strategies are significantly different.We also calculated and discussed the differences between our modeled height map and the two reference products in a view of the different canopy height modeling approaches.

Canopy Height Map in China
A contiguous map of canopy heights with a spatial resolution of 1 km over the continental China was generated from the FFNN using the gridded geospatial predictors (Figure 3a).This map (mean = 36.44m, Figure 3a) showed a good consistency with the spatial pattern of GLAS tree heights (mean = 30.00m, Figure 3b).From the estimation result of canopy heights over China, we can see that relatively tall trees were growing in the central and southern regions (Sichuan, Hubei, Yunnan and Chongqing).The trees distributed in the Northern China show obviously lower heights than those in Southern regions.Specially, the forests distributed in Heilongjiang and East of Inner Mongolia province, which are close to the border of Russia, have the lowest tree heights.The varying distribution of tree heights over China likely captures climate gradients.However, further assessment is still needed for exploring the underlying patterns of the tree height distribution.

Ground Validation and Error Analysis
92 field plots were used in comparison with the FFNN-modeled tree heights.The validation showed a good agreement between the modeled and field-measured tree heights (Figure 4a).In order to show the outperformance of our model with respect to previous studies, we also evaluated HSimrad, HNi and the GLAS tree height map produced from 20 GLAS shots selected according to Equation (2).From these comparison results (Figure 4b-d), we can see that the modeled heights have better accuracy than other approaches.

Actual GLAS-Derived Tree Height Validation and Error Analysis
The modeled heights from the FFNN at the continental China scale were also evaluated using the GLAS-derived heights.This validation was based on two types of GLAS metrics, including the footprint-level GLAS heights and a full gridded GLAS height map.Both types of validation showed a consistency agreement between the FFNN modeled and GLAS height metrics (R = 0.725, RMSE = 7.86 m for the footprint-to-pixel level; R = 0.759, RMSE = 8.85 for the pixel-to-pixel level; Figure 5).However, Figure 5a depicts some difference between the model tree heights and footprint-to-pixel-level tree heights derived from GLAS shots.The main reason is that tree heights model in this study were trained from 20 GLAS shots which were selected according to Equation (2), including 15 training data and five validation data.As for one GLAS shot location ( ), these 20 GLAS shots for training and validating the tree height model should have the most similar temperature and precipitation with the GLAS shot location ( ).Additionally, the GLAS shot of location was included in the 20 GLAS shots.In general, modeled tree height should be similar to the GLAS shot height of location .When the temperature and precipitation of the GLAS shot location ( ) have a significant difference from the other 19 GLAS shots, the corresponding modeled tree height might be different from the GLAS shot height of location .Conversely, the modeled tree height is still similar to the mean value height of these 20 GLAS shots, which can explain the better consistency agreement in Figure 5b than in Figure 5a.
In order to further evaluate the modeled tree heights, we demonstrated the spatial distribution of deviations between the modeled and GLAS-derived heights.As shown in the map depicting the difference at the pixel-to-pixel level (Figure 6a), most forested pixels showed a moderate predictive power but overestimations were dominant over some regions.The overestimation of tree heights tended to distribute over some special area where temperature is relatively low and precipitation is relatively high (e.g., east of Heilongjiang and Jilin, southeast of Tibet).The overestimation may be caused by the tree height model, which is more sensitive to the precipitation than temperature during the model training process.In addition, relatively high precipitation may have bigger impacts on forest trees' growing and characteristics' expression (leaf reflectance of forest tree) than a relatively low temperature, which might affect the accuracy of modeled tree heights.The specific reasons will be explored in a further study.

Comparison with Simard Tree Heights Map
Based on the definition of forest land area in this study, we obtained the pixel-to-pixel comparison result between our modeled heights and the .The difference map (Figure 7a) showed a moderate agreement.The differences were less generally 10 m.
From the comparison result (R = 0.738, RMSE = 7.65 m; Figure 7b), we can see a very good consistency between modeled tree heights and the when tree heights are less than 40 m.It is noteworthy that there are some obvious discrepancies existing where modeled tree heights are higher than 40 m.The disagreement between modeled tree heights and the for trees greater than 40 m might be caused by the following two reasons.On the one hand, the Simard tree heights are is the regional average height which made the forest height prediction result likely being less than 40 m.On the other hand, the overestimation of the modeled tree heights above 40min some places also appears to be an obvious disagreement with the .

Inter-Comparison with Ni Tree Heights
Based on the climate zones defined by Ni [3], the comparison between the modeled heights and the HNi was conducted (Figure 8).As depicted in Figure 8a, the modeled heights and the HNi generally agreed.For the most of forest land area, the difference between the modeled heights and the HNi was less than 10 m.However, some regions in Southern China had relatively more significant difference, which might be caused by higher HNi values in this region.The same comparison result was shown in Figure 8b (R 2 = 0.519, and RMSE = 8.99 m).In the case that the heights are less than 40 m, there is a good agreement between the modeled heights and the HNi.However, when the heights are higher than 40 m, the difference became larger.The main reason for these obvious errors is that in this process of comparison, we calculated the mean value of the modeled heights in climate zones, while the HNi is the maximum tree heights predicted from the optimized allometric scaling and resource limitation (ASRL) model as the height of the climate zone.

Concluding Remarks
The objective of this work was to predict the large-scale spatial pattern of forest heights over the continental China.The GLAS derived heights were used as an input for the model and as a result, a spatially-continuous forest height map was produced with a spatial resolution of 1000 m.Additional geospatial predictors were climatic variables (temperature and precipitation), forest type, tree cover percent and land surface reflectance.The approach was designed to train the artificial neural network (ANN) model to minimize the errors between GLAS-derived actual heights and model-derived forest heights.For each pixel of the forested land defined in the study, there was one ANN model to be trained and tested.
The ANN-modeled heights were firstly evaluated at the site scale.In three field experimental sites, regression analysis was performed to assess the correlation between the modeled and field-measured tree heights.The comparison showed a good agreement (n = 92, R = 0.827, and RMSE = 4.15 m).Meanwhile the modeled heights were assessed at the GLAS footprint-to-pixel level.We used more than thirty thousand GLAS shots for this evaluation.The comparison between the modeled and GLAS heights showed a reasonable correspondence (R = 0.725, and RMSE = 7.86 m).
In addition, the ANN modeled heights were compared with continuous forest height maps over the China continent.The comparison was firstly conducted with the interpolated GLAS height map.Their difference map showed that there were relatively small errors in the most Chinese forested lands (<10 m), and the ANN model tended to predict little higher tree heights than the gridded GLAS heights.Furthermore, the modeled tree heights were validated using the published products of tree heights (RMSE = 7.65 m, R = 0.738 for HSimard; R 2 = 0.519, RMSE = 8.99 for HNi).According to the difference maps between the modeled and existing products, we obtained a small overestimation compared with the HSimard and an underestimation compared with the HNi.These discrepancies could be mainly attributed to some inherent limitations of various tree height estimation approaches.
This paper reported the ANN model for the extensive forest canopy height predictions.The training method of the ANN model for each pixel successfully compensated for the limitations of previous tree height estimation approaches by using a single model in large-scale regions, which can greatly improve the accuracy of forest tree height estimation to some extent.As shown in the results, the trained ANN model in this article demonstrated a good performance in the estimation of tree heights over the continental China.
Nevertheless, the ANN model is a black-box model, which cannot describe the growth mechanism of trees.For this reason, the trained ANN model in this paper still yielded ambiguous results in some regions.These ambiguous results were possibly caused by the uncertainties of input data in the areas with complex geographic and environmental conditions, or by the inaccurate training data due to the topographic effects.Further studies should focus on reducing the impacts of environment and topography on the uncertainties of tree height estimation, and extending the approach to biomass prediction.

Figure 1 .
Figure 1.Field measurement sites.(a) Distribution of three field measurement regions; (b) Dayekou field sites in Gansu; (c) Taihe field sites in Jiangxi; (d) Puer field sites in Yunnan; (b-d) the background color map represents the MODIS land cover (LC) forest types.
depicts the schematic diagram of the ANN model proposed in this study.

Figure 2 .
Figure 2. The schematic diagram of the ANN model proposed in this study.There are 11 neurons in the input layer, 11 neurons in the hidden layer and 1 neuron in the output layer.

Figure 3 .
Figure 3. Maps of modeled and averaged GLAS-derived tree heights over continental China at 1-km spatial resolution.(a) The modeled height map using the trained ANN model proposed in this paper; (b) GLAS tree height map.The GLAS tree height map is made from the interpolation of the averaged GLAS-derived tree heights.

Figure 4 .
Figure 4. Ground validation result of the modeled tree heights, HSimrad, HNi and GLAS tree height map.(a) The comparison between modeled tree heights and field measured tree heights; (b) The comparison between HSimrad and field measured tree heights; (c) The comparison between HNi and field measured tree heights; (d) The comparison between averaged GLAS tree heights and field measured tree heights.

Figure 5 .
Figure5.The comparison results between modeled tree heights and actual GLAS derived tree heights.(a) The comparison between modeled tree heights and discrete GLAS tree heights; (b) the comparison between modeled tree heights and averaged GLAS tree heights; (c) statistical characteristics of modeled tree heights, discrete GLAS tree heights and averaged GLAS tree heights.

Figure 6 .
Figure 6.The difference map of modeled tree heights with averaged GLAS tree heights.

Figure 7 .
Figure 7.The comparison result between modeled tree heights and Simard tree heights.(a) The difference map of modeled tree heights with Simard tree heights; (b) the comparison between modeled tree heights and Simard tree heights.

Figure 8 .
Figure 8.The comparison result between the modeled tree heights and the HNi s.(a) The difference map of the modeled tree heights against HNi; (b) the comparison between the modeled tree heights and HNi.

Table 1 .
The detailed information on three field experiment regions.