ESTIMATION OF REGIONAL FOREST ABOVEGROUND BIOMASS COMBINING ICESAT-GLAS WAVEFORMS AND HJ-1 A / HSI HYPERSPECTRAL IMAGERIES

Estimation of forest aboveground biomass (AGB) is a critical challenge for understanding the global carbon cycle because it dominates the dynamics of the terrestrial carbon cycle. Light Detection and Ranging (LiDAR) system has a unique capability for estimating accurately forest canopy height, which has a direct relationship and can provide better understanding to the forest AGB. The Geoscience Laser Altimeter System (GLAS) onboard the Ice, Cloud, and land Elevation Satellite (ICESat) is the first polarorbiting LiDAR instrument for global observations of Earth, and it has been widely used for extracting forest AGB with footprints of nominally 70m in diameter on the earth’s surface. However, the GLAS footprints are discrete geographically, and thus it has been restricted to produce the regional full coverage of forest AGB. To overcome the limit of discontinuity, the Hyper Spectral Imager (HSI) of HJ-1A with 115 bands was combined with GLAS waveforms to predict the regional forest AGB in the study. Corresponding with the field investigation in Wangqing of Changbai Mountain, China, the GLAS waveform metrics were derived and employed to establish the AGB model, which was used further for estimating the AGB within GLAS footprints. For HSI imagery, the Minimum Noise Fraction (MNF) method was used to decrease noise and reduce the dimensionality of spectral bands, and consequently the first three of MNF were able to offer almost 98% spectral information and qualified to regress with the GLAS estimated AGB. Afterwards, the support vector regression (SVR) method was employed in the study to establish the relationship between GLAS estimated AGB and three of HSI MNF (i.e. MNF1, MNF2 and MNF3), and accordingly the full covered regional forest AGB map was produced. The results showed that the adj.R and RMSE of SVR-AGB models were 0.75 and 4.68 t·hm for broadleaf forests, 0.73 and 5.39 t·hm for coniferous forests and 0.71 and 6.15 t·hm for mixed forests respectively. The full covered regional forest AGB map of the study area had 0.62 of accuracy and 11.11 t·hm of RMSE. The study demonstrated that it holds great potential to achieve the full covered regional forest AGB distribution with higher accuracy by combing LiDAR data and hyperspectral imageries. * Yanqiu Xing: the corresponding author.


INTRODUCTION
Estimating and preserving the carbon stock in forests can help devise sequestration strategies and reduce greenhouse gas emissions (Zhang, 2014;Canadell, 2007;Miles, 2008).Forest aboveground biomass (AGB) has received more and more attention during the last decades due to its relevance to carbon cycle balancing and species diversity (Saatchi,2011a;Laurin, 2014;Saatchi, 2011b).Therefore, accurate calculation of Forest AGB and their distribution is of great significance to research on carbon cycles and carbon stocks (Neigh, 2013).
Light Detection and Ranging (LiDAR) has a unique capability for estimating forest structure parameters, especially forest vertical structure parameters.Airborne small-footprint LiDAR is considered the most accurate remote sensing technology for mapping forest AGB (Laurin, 2014;Zolkos, 2013), but its disadvantages such as expensive costs and complex data processing, limit its application on large-scale areas.On the other hand, space-borne LiDAR has a wide observation scope and can capture large-scale even global information without the limitation of time and weather (Asner, 2010;Drake, 2002).
The Geoscience Laser Altimeter System (GLAS) onboard the Ice, Cloud, and land Elevation satellite (ICESat) launched on 12 January 2003, a unique LiDAR instrument designed for continuous global observation of the earth (Zwally, 2002;Sun, 2008).It is able to capture full waveform data consists of the first return waveform reflected by forests and the last return waveform reflected by the ground.Therefore, we can obtain vertical distribution information such as forest height through analysing the full waveform data.The GLAS waveform has been widely used to estimate forest structure parameters and forest AGB (Rosette, 2008;Pflugmacher, 2008;Xing, 2010;Lefsky, 2005;Enßle, 2014;Sun, 2008).Although the GLAS performs well in estimating forest structure parameters, but its sparse distribution is discontinuity which limits its application on large-scale areas alone.The combination of multiple sources of data is required to map large areas of forest AGB (Sawada, 2015;Chi, 2015;Guo, 2010).In addition, GLAS waveform is sensitive to terrain slope due to its large footprint size, which causes overestimate of forest height and decreases forest AGB estimation accuracy.
In this paper, in order to overcome the two problems, ICESat-GLAS waveform metrics (W, I, TS) were combined with three MNF components of HJ-1A/HSI hyperspectral imageries to estimate forest AGB based on Support Vector Regression (SVR) method in order to improve the forest AGB estimation accuracy and achieve the full covered regional forest AGB distribution.

Study area
The study area is within Wangqing forestry of Jilin province, China, and is located along the border between China and North Korea (43°05′N-43°40′N, 129°56′E-131°04′E) (see Figure 1).It belongs to Changbai mountain system, one of the most valuable reserves in China, and it is a rich gene pool of many plant species.The elevation ranges from 360m to 1477m, and the terrain slope generally ranges from 0° to 45°.

ICESat-GLAS waveform
ICESat-GLAS carries three lasers in total, and each laser emits 40 pulses of 1064 nm per second, resulting in an elliptical footprint of ~70m diameter on the surface of the earth, with each footprint separated by 172m (Schutz, 2005).GLAS totally offers fifteen kinds of products (Sun, 2008).In the study, release-31 of the GLA01 and GLA14 data were acquired through the National Snow and Ice Data Center (NSIDC) website (http://nsidc.org/data/icesat/).The GLA01 data contains raw waveform data.For land surfaces, the waveform has 544 bins with one bin represents 15cm.The GLA14 data does not contain waveform data, but includes information about the observation conditions and the latitude and longitude of footprints, etc.The two products are joined together by the record index of the GLAS footprints.
Since GLAS waveforms was easily affected by clouds and system noise, cloud-contaminated and abnormal waveforms were removed before extracting waveform parameters by using the saturation-free flag (i_satNdx=0) and cloud-free flag (i_FRir_qaFlg=15) data extracted from GLA14 product (Popescu, 2011).

HJ-1A HSI imageries
The Disaster and Environment Monitoring and Forecast constellation consists of two satellites (HJ-1A and HJ-1B), and was launched on September 6, 2008.HJ-1A is equipped with a Hyper Spectral Imager (HSI) with 50km swath and 100 m spatial resolution consisting of 115 bands with 4.32nm spectral resolution in the 0.45~ 0.95μm range.Four images of level-2 data were acquired through the China Centre for Resources Satellite Data and Application web site and used for Forest AGB estimation (http://www.cresda.com/n16/index.html).

Field inventory data
The field data collection was carried out in September of 2006 and 2007 in the study area and 183 circular plots of 500 m 2 along ICESat-GLAS flight track direction were randomly established.The description of field data was listed in Table 1.These plots consist of 101 broad-leaved forest, 38 coniferous forest and 44 mixed forest.In each plot, we recorded species and forest canopy density.Tree height and DBH (Diameter at Breast Height) were also measured with a DBH greater than 10 cm.For trees with a DBH between 4cm and 10cm, we only recorded tree species, and we made no record of trees with a DBH less than 4 cm.In addition, the terrain slope and the number of trees were also recorded in each plot.Biomass of single trees in each plot was calculated using the polynomial function with DBH as the variable developed by Deo (Deo, 2008).In the study, the sum of the single tree biomass divided by plot area was defined as the forest aboveground biomass of each plot.

METHOD
3.1 Data processing

ICESat-GLAS waveform processing and metrics extraction:
The GLAS waveforms were smoothed by Gaussian filter to remove noise.The smoothed waveform was modelled with Gaussian components using the algorithm developed by Brenner et al (Brenner, 2003).The signal start and signal end were identified using the background noise threshold which was set to the mean background noise plus 4.5 times the background noise standard deviation (Lefsky, 2005).
Generally, waveform extent (W) as shown in Figure 2 was defined as maximum forest height on a flat area.But the waveform was increased as the function of terrain slope and footprint size in mountain forests.To increase the accuracy of maximum forest height and forest AGB, the terrain slope parameter TS was extracted from the GLAS waveform based on an approach developed by Mahoney (Mahoney, 2014) Where TS =terrain slope W gf =the of the Gaussian waveform corresponding to ground return W m =the minimum width of the Gaussian waveform caused by the duration of the emitted signal and atmospheric attenuation D =the diameter of the GLAS footprint A =the GLAS waveform maximum amplitude.
In this paper the forest canopy density was calculated for each footprint as the ratio of the forest canopy energy to the total energy of the GLAS waveform and the equation is shown in Equation (2).As shown in figure 2, the forest canopy energy was calculated by summing the reflected energy of forest canopy height bins of the GLAS waveform, and the total energy was computed by adding up the energy of the GLAS waveform from signal start to signal end.
Where I = forest canopy density E v =the reflected energy of forest canopy height bins of the GLAS waveform E =the total energy of GLAS waveform.Because HSI hyperspectral image consists of 115 bands, there is strong data redundancy and high correlation between neighboring bands.In the study, the Minimum Noise Fraction (MNF) method was applied to decrease noise and reduce the dimensionality of spectral bands.The first three MNF components (i.e.MNF1, MNF2 and MNF3) contains almost 98% of total information and satisfies research requirement and so they were used to do forest AGB estimation.Considering that the spatial resolution of GLAS data and HSI data are different, circle buffers of 70m in diameter were established on the HSI image with the coordinate of GLAS footprints as their center, and the mean value of pixels covered by circle buffers was defined as MNF value of the corresponding GLAS footprints.

Regional forest AGB estimation combining GLAS waveform and HJ-1A/HSI imageries
In the study, GLAS waveform and HSI image were combined to estimate forest AGB.Metrics (W, TS, I) derived from GLAS data and field AGB were used to develop the GLAS-AGB models by multiple regression method for different forest types in the study area.The GLAS-AGB models were then used to calculate forest AGB of the remaining GLAS footprints.
The SVR (Support Vector Regression) method was applied to combine the GLAS metrics with HSI metrics to finish regional forest AGB estimation.SVR is based on statistical learning theory and the structural risk minimization principle.It performs well in approximating and generalizing.The basic principle of SVR is to map the data into a high dimensional feature space via a kernel function, after which a linear regression is performed in this feature space.The radial basis function is a normal kernel function as shown in Equation(3): 2 ( , ) exp( ) Where  = the width of the kernel function x i =the centre of the kernel function x =the dependent variable In the paper, we used the LIBSVM tools for the forest AGB estimation.To build the SVR-AGB model, all the forest AGB of GLAS footprints and the correlated three MNF components of HSI data were put into LIBSVM, and they were randomly divided into training and validation data to develop SVR-AGB models for three forest types.The SVR-AGB models were then used to calculate forest AGB of each pixel of HSI data and the AGB map of the entire study area was achieved.The flow chart of the forest AGB map process is shown in Figure 3.
The adjusted coefficient of determination (adj.R 2 ) was used to evaluate the goodness of fit for the model.The root mean square error (RMSE) was used to assess the error between the estimated AGB and field-observed AGB at the footprint level.

Biomass estimation from GLAS waveform parameters
As shown in Figure 4, TS has a strong linear relationship with terrain slope (adj.R 2 =0.78,RMSE=4.74°).Therefore, TS can be used to correct the influence of terrain slope on GLAS waveform for further forest AGB estimation.
As shown in table 2, when the under-story vegetation height was set as 2m, the accuracy was the highest (adj.R 2 =0.64,RMSE=0.13).Thus in this study, the GLAS parameter I extracted from the GLAS waveform when the under-story vegetation height was set as 2m was used to estimate forest AGB.Additionally, from the Table 3, it can also be seen that the adj.R 2 value of AGB models with W and I as independent variables of the three forest types significantly increased, demonstrating that forest AGB really has linear relationship with I.When all of the GLAS parameters (W, TS, I) were considered in the AGB models, the accuracy was the best.In this case, adj.R 2 and RMSE values of the forest AGB model are 0.65 and 10.21 t•hm -2 for broadleaf forests, 0.66 and 11.86 t•hm -2 for coniferous forests and 0.62 and 12.85 t•hm -2 for mixed forests.
The regression equations and the validation results of AGB models of the three forest types are shown in Table 4.The results showed that the forest AGB estimated from GLAS waveforms and the forest AGB calculated from forest inventory data has good consistency with adj.R 2 =0.75 for broadleaf forest, adj.R 2 =0.77 for coniferous forest and adj.R 2 =0.76 for mixed forest.Therefore, GLAS-AGB estimation models can be used to predict regional forest aboveground biomass combining with HJ-1A/HSI data in further study.

Regional AGB estimates and validation
The regional SVR-AGB models are shown in Table 5.It can be seen that the accuracy of broadleaf forest was higher than coniferous forest and mixed forest.The forest AGB of the study area was calculated using the SVR-AGB model pixel-by-pixel.
The final results are shown in Figure 5, and it could be seen that the forest AGB of the study area ranges from 0 to 164 t•hm -2 .The highest forest AGB is located in the northern and southwestern part of the study area where evergreen trees are dominant.The southern and northeastern parts of the study area are dominant with larch whose leaves had become yellow or fallen, resulting in a lower AGB value.The regions of the study area with a forest AGB value of zero are mainly residential areas, roads and rivers.
Even though the accuracy of forest AGB estimation was improved, there is also some bias between the predicted forest AGB and field forest AGB, and the reasons are various.As shown in Figure 5, the forest AGB increases as altitude increases and then declined when the altitude reached 700m.Large amounts of forest AGB were concentrated in the altitudes between 400m and 800m.Because the study area is dominated by a cool temperate continental climate, when the altitude reached 1000m, the temperature is only 3 , which is not ℃ suitable for plants, and there were only some small but old trees.The time of acquisition of the data in the study was different, which may also bring some negative influence on the results.Table .5Results of forest AGB models estimated with HJ-1A/HSI parameters 0 -1 0 0 1 0 0 -2 0 0 2 0 0 -3 0 0 3 0 0 -4 0 0 4 0 0 -5 0 0 5 0 0 -6 0 0 6 0 0 -7 0 0 7 0 0 -8 0 0 8 0 0 -9 0 0 9 0 0 -1 0 0 0 1 0 0 0 -1 1 0 0 1 1 0 0 -1 2 0 0 The results demonstrated that GLAS waveform plays an important role in linking field inventory data to HJ-1A/HSI data, because of the fact that the acquisition of field inventory data is both time-and resource consuming.GLAS waveform records the interaction of laser beams with forest and ground which reflect the vertical structure of forest.Some GLAS parameters such as W, I, and TS can be used to develop a forest AGB model.TS has a strong relationship with terrain slope, and can be used to correct the influence of terrain slope on the GLAS waveform.The ratio I of the forest canopy energy to the total energy of GLAS waveform has the best correlation with forest canopy density when the under-story vegetation height was 2m.And the results showed that GLAS waveform parameters have the good relationship with AGB and the AGB equations for the three forest types can be used to obtain sample biomass data for regional AGB estimation using HJ-1A/HSI images.
GLAS parameters which reflecting forest vertical and horizontal structure are considered in the AGB model, and the accuracy of forest AGB model is improved, comparing with that alone using vertical structure or horizontal structure parameters.
From the results shown in Table 5, we can see that the combination of GLAS data and HJ-1A/HSI data is reasonable and promising for regional forest AGB estimation.Combining GLAS and HJ-1A/HSI data can improve the estimation accuracy of forest AGB and can complete regional forest AGB estimation.However, because of the spatial resolution of HJ-1A/HSI data, the spatial distribution and the limitation number of GLAS footprints, the accuracy of the estimated forest AGB of the study area was not very high.In the future studies, some high spatial resolution and high spectral resolution image data can be considered to improve the accuracy.

Figure 1 .
Figure 1.Study area and the location of sample points Figure 2. The GLAS waveform and GLAS parameters

Figure 3 .
Figure 3. Flow chart of regional biomass estimation from GLAS, HJ-1A/HSI, and field data 4 RESULTS AND DISCUSSION

Figure 4 .
Figure 4. Scatter plot of TS and terrain slope

Figure 4 .
Figure 4. Map of Forest aboveground biomass Figure 5. Forest aboveground biomass with different altitude

Table . 1
Description of field data

Table . 2
Adj.R 2 and RMSE value of GLAS forest canopy density model with different under story vegetation height

Table . 4
Regression equations of AGB models of the three forest typesAs shown in Table3, when W and TS were considered in the regression analysis, adj.R 2 significantly increased and RMSE significantly declined.For example, adj.R 2 increased from 0.39 to 0.45 in broadleaf forests, from 0.36 to 0.53 in coniferous forests and from 0.46 to 0.52 in mixed forests.RMSE declined from 18.64 t•hm -2 to 16.23 t•hm -2 in broadleaf forests, from 20.25 t•hm -2 to 16.61 t•hm -2 in coniferous forests, and from 19.56 t•hm -2 to 15.19 t•hm -2 in mixed forests, respectively.