National Forest Aboveground Biomass Mapping from ICESat/GLAS Data and MODIS Imagery in China

Abstract

Forest aboveground biomass (AGB) was mapped throughout China using large footprint LiDAR waveform data from the Geoscience Laser Altimeter System (GLAS) onboard NASA’s Ice, Cloud, and land Elevation Satellite (ICESat), Moderate Resolution Imaging Spectro-radiometer (MODIS) imagery and forest inventory data. The entire land of China was divided into seven zones according to the geographic characteristics of the forests. The forest AGB prediction models were separately developed for different forest types in each of the seven forest zones at GLAS footprint level from GLAS waveform parameters and biomass derived from height and diameter at breast height (DBH) field observation. Some waveform parameters used in the prediction models were able to reduce the effects of slope on biomass estimation. The models of GLAS-based biomass estimates were developed by using GLAS footprints with slopes less than 20° and slopes ≥ 20°, respectively. Then, all GLAS footprint biomass and MODIS data were used to establish Random Forest regression models for extrapolating footprint AGB to a nationwide scale. The total amount of estimated AGB in Chinese forests around 2006 was about 12,622 Mt vs. 12,617 Mt derived from the seventh national forest resource inventory data. Nearly half of all provinces showed a relative error (%) of less than 20%, and 80% of total provinces had relative errors less than 50%.

1. Introduction

As the major component of terrestrial ecosystems, forests hold 70%–80% of terrestrial aboveground and belowground biomass and play an essential role in the global carbon cycle and climate change [1]. Forests are vulnerable to fire, logging, pests, and land conversion, which release carbon easily to the atmosphere [2]. The annual carbon flux between forests and the atmosphere accounts for 90% of the flux between the atmosphere and Earth’s land surface [3]. Therefore, mapping the magnitude, spatial distribution, and change of forest aboveground biomass (AGB) over time is important for improving estimates of terrestrial carbon sources and sinks. Furthermore, some studies of global change and the carbon cycle showed that there was a large uncertainty in the forest carbon stocks in tropical and boreal regions because of the lack of accurate forest biomass maps [4,5]. Regional to global information on the human impact on carbon stocks and ecological balance also requires accurate determination of forest biomass and monitoring of its changes, especially in areas of fragmented forest cover and in developing countries [6,7,8]. These scientific objectives require accurate mapping of forest biomass at the national level for the assessment of forest carbon stocks and their dynamics.

Traditionally, field measurements for estimation AGB are collected at sampling sites by measuring tree heights and DBH (diameter at breast height). The biomass of a tree can be estimated from its specific allometric equation using these two variables or DBH alone. The biomass of a plot is typically expressed in the form of biomass density, i.e., in weight per unit area. Previous studies on biomass estimation were developed based on the extrapolation of biomass from specific sites to large areas (i.e., multiplying the mean biomass density estimated from direct field measurements by the forested area or using biomass expansion factor to convert timber volume to biomass) [2,9,10,11]. These methods can obtain good biomass estimation for a small area, such as a forest stand. However, the lack of field data in remote areas and the inconsistency of data collection methods among different management units over large regions were the major constraints to obtaining reliable large-scale biomass estimation using field-based methods [12,13]. Moreover, obtaining comprehensive, spatially complete, temporally uniform, and accurate forest inventory data was usually time-consuming and labor-intensive over huge areas and field campaigns were not well suited for detecting changes because a single measurement campaign can extend over several years in most countries, especially in developing countries with large land area [12,14,15].

Satellite remote sensing is fairly advantageous for regional or global AGB estimations because of its large area coverage [16]. However, as no remote sensing instrument has been developed that is capable of providing direct measurement of biomass, additional field-sampled biomass is required for establishing relationships between remote sensing signals and biomass [17]. Optical satellite data can provide fine results in the discrimination of largely different biomass classes across a variety of spatial and temporal scales due to its spectral sensitivity to different species. Furthermore, satellite-based optical imaging spectrometers are perceptive to vegetation structural parameters such as tree abundance, basal area, stem density and crown size, which are correlated to some degree with biomass because vegetation spectral reflectance contains information on the vegetation chlorophyll absorption bands in the visible region (especially in blue band and red band) and sustained high reflectance in the near-infrared region [18]. Particularly, recent research has focused on the utility of short wave infrared (SWIR) data for estimating AGB [19,20]. Besides application of single-band spectral signatures, various vegetation indices (VIs) derived from TM (Thematic Mapper), ETM+ (Enhanced Thematic Mapper Plus), AVHRR (Advanced Very High Resolution Radiometer), and MODIS data have been used to estimate forest biomass in various areas, including southern Asia, boreal and temperate forest zones in America and Asia, subtropical areas in South America, southwestern United States and tropical Africa, America and Asia [4,19,20,21,22,23,24,25]. Although passive optical remote sensing is a less expensive technology for biomass estimation over vast areas, frequent cloud coverage in mountain areas and the low saturation level of spectral bands and VIs on biomass estimation in medium to high biomass forests are the major disadvantages of this technology [26,27,28].

LiDAR (light detection and ranging) is an active sensor type based on laser ranging, which emits a laser pulse and measures the distance based on half the elapsed time between the emission of a pulse and the detection of a reflected return [29]. When the surface is covered by vegetation, echoes are a function of the vertical distribution of vegetation and the ground surface within a footprint. LiDAR data capture the horizontal and vertical structure of vegetation in great detail, resulting in accurate estimates of forest biomass across a broad range of forest types and biomes [30,31,32,33]. LiDAR systems acquire data either over small footprints (<1 m, point-cloud data, or waveform data, airborne) or large footprints (>10 m, waveform, either airborne or spaceborne) [34,35]. LiDAR is the only sensor type available at present whose signal does not saturate in high biomass forests (e.g., 1200 Mg·ha⁻¹ and 1300 Mg·ha⁻¹ [36,37]). Estimating biomass with airborne LiDAR data is often more accurate [38,39], but the associated large data volume, as well as sophisticated technical equipment and high acquisition costs to observe remote areas limit its application at regional and global scales [16,19,27]. The Geoscience Laser Altimeter System (GLAS), onboard the Ice, Cloud, and land Elevation Satellite (ICESat), was a spaceborne LiDAR system recording full waveforms over large footprints [40]. GLAS used 1064-nm laser pulses operating at 40 Hz and recorded the echo of those pulses from footprints ~65 m in diameter, spaced 172 m apart [41]. GLAS footprint data as a manner of sampling are crucial for associating field plot data with optical imaging systems because LiDAR samples could compensate for the lack of systematic spatial sampling of AGB from ground measurements. One major limitation of GLAS was the lack of imaging capability and the fact that it provided relative sparse sampling information on forest structure [42]. Therefore, multi-sensor data synergy is required to estimate forest structural parameters at the regional scale.

Over the past few years, the integration of optical remote sensing and GLAS data has been successfully used to estimate forest volume and biomass because the signature from these two types of data provide relevant information content, e.g., species and forest coverage from the optical data, and three-dimensional vertical forest structure information from LiDAR data [42,43]. For example, Boudreau et al. [12] used GLAS waveform metrics, airborne profiling LiDAR, and field plots to estimate forest volume, biomass, and carbon in Quebec, Canada, covering an area of 1.3 million km². A generic airborne LiDAR-based biomass equation (R² = 0.65) was developed from field inventory data as a function of airborne profiling metrics. A second model was explored that predicted biomass as a function of GLAS waveform variables. Thus, GLAS footprints became the regional sampling tool used to estimate forest biomass on forest cover strata that were derived from Landsat-7 ETM+ [44]. Nelson et al. [30] developed a GLAS-based equation to estimate timber volume in south-central Siberia using GLAS and MODIS data. Although only 51 ground plots co-located with GLAS footprints were sampled, regional total volume estimates and per-hectare estimates were compared with ground-based study results for the entire 811,414 km² area, with differences of 1.1% and 11.9%, respectively. Baccini et al. [45] mapped AGB across tropical Africa (1-km resolution) using MODIS data with a large dataset of field measurements. They found a strong relationship between the estimated AGB and GLAS height metrics (average height and height of median energy). That study noted, “…there are currently limited high quality field biomass estimates available at sufficient spatial extent to develop and independently validate maps of AGB across tropical regions…” [45]. Baccini et al. [20] estimated the carbon density of aboveground live wood vegetation for the pan-tropics at a spatial resolution of 500 m from MODIS data using GLAS waveform metrics. Saatchi et al. [4] mapped biomass (above- and belowground) in tropical regions across three continents using a combination of data from 4079 in situ inventory plots and GLAS samples of forest structure, plus optical and microwave imagery (1-km resolution). They developed continent-based allometric equations to provide the best models to convert Lorey’s heights derived from GLAS data to forest biomass. A data fusion model based on the maximum entropy approach was used to extrapolate AGB derived from GLAS footprints to the entire landscape. This benchmark map provided comparable estimates of AGB for 75 developing countries.

China, as one of the world’s fastest developing countries, needs to produce robust estimates of forest biomass and carbon stocks for successful implementation of climate change mitigation policies. As one of the five most forest-rich countries [46], China is rich in temperate forests and subtropical forests. Timely and accurate measurements of forest biomass and its distribution are increasingly needed to support a wide range of activities related to sustainable forest management and carbon accounting. Previous studies on estimates of forest biomass in China were based on statistical analysis of the biomass-volume relationship based on nationwide forest inventory data [47,48]. Despite the high precision of such inventories, they do not provide maps of biomass at a resolution useful for assessing land-use change. An AGB map of China with clear and detailed spatial distribution is urgently needed. However, the benchmark map of Saatchi et al. [4] did not cover the entire land of China. Similarly, the pan-tropical map generated by Baccini et al. [20] only covered the area in southern China (below 30°N). More importantly, no field survey samples from the Chinese territory were included in the two studies.

In this study, we explore the capabilities and limitations of satellite remote sensing data (GLAS and MODIS) and ground measurements for mapping AGB in China. The slope effect on biomass estimation in a footprint was considered. The potential information on biomass from GLAS waveforms and GLAS samples for reliable biomass estimation were studied using field data. Proper prediction models were used to extend the biomass estimates from GLAS samples to all forested areas throughout China. Biomass estimation results were compared with field biomass data within a footprint, with the aim of analyzing the effect of slope on the accuracy of the estimates. The biomass map was validated with other biomass estimates derived from forest inventory data.

2. Materials and Methods

2.1. Study Area

The territory of China lies approximately between 18°N and 54°N, and 73°E and 135°E, with a total area of approximately 9,600,000 km². China is a mountainous country, where plateaus and mountainous regions account for about 60% of the total land area. Plains account for nearly one-fifth, and the rest consists of basins and hills [49]. China’s climate is dominated by a continental monsoon climate. Seasonal changes and annual variability of temperature and precipitation are significant in most regions of China, which are major factors in the formation of complex and diverse climate as well as topography.

Figure 1. Seven forest zones in China: Red boxes are the field sample sites. (I = cold temperate zone; II = temperate zone; III = warm temperate zone; IV = subtropical zone; V = tropical zone; VI = Neimeng-Xinjiang arid zone; and VII = Qinghai-Xizang plateau alpine zone).

Figure 1. Seven forest zones in China: Red boxes are the field sample sites. (I = cold temperate zone; II = temperate zone; III = warm temperate zone; IV = subtropical zone; V = tropical zone; VI = Neimeng-Xinjiang arid zone; and VII = Qinghai-Xizang plateau alpine zone).

Precipitation decreases from the coast to the inland areas. Three vegetation regions, namely forests, steppes, and deserts, correspond to moist, semi-arid, and arid climates, respectively [50]. According to Chinese vegetation geography [49,51] and geographic characteristics of forests, Chinese forests have been divided into five forest zones (Figure 1): Cold temperate zone (I), which is dominated by deciduous needleleaf forests; temperate zone (II), which is characterized by deciduous mixed broadleaf-needleleaf forests; warm temperate zone (III), including China’s largest plain, with secondary broadleaved-mixed forests and intensive agricultural activity; subtropical zone (IV), with large formations of evergreen broadleaf forests (the western region of the subtropical zone is dominated by high mountains and affected by the southwest monsoon, while the eastern region is dominated by hills and affected by the southeast monsoon); tropical zone (V), whose annual average temperature is over 22 °C and average annual precipitation is above 1500 mm. Two other vegetation zones (not dominated by forests) were also used in this study (Figure 1): The Neimeng-Xinjiang arid zone (VI), which is distinguished by Picea and Larix in the Tianshan, Altai, and Qilian Mountain regions (most land areas in this zone are covered by steppes and deserts due to a continental climate with severe annual variations in temperature); and the Qinghai-Xizang plateau alpine zone (VII), which retains the largest area of virgin forest in China.

2.2. Data

2.2.1. Field Data

Ground survey data spanning a variety of forest types in each forest zone (except the Neimeng-Xinjiang arid zone) were obtained from two sources, as listed in Table 1. A database of forest inventory data was assembled from the National Forest Resource Inventory (NFRI) of the State Forestry Administration, China. The survey results of NFRI data were shown in the form of forest maps, which were drawn from the statistical stratification procedure of forest stands. A stand is a contiguous area that contains a community of trees that are relatively homogeneous or have a common set of characteristics (e.g., age and size class distribution, composition, structure, spatial arrangement, site quality, condition, or location to distinguish it from adjacent communities) [52]. Normally, a stand is studied or managed as a single unit. The field investigation includes records of stand area, forest types, dominant species, stem volume, stand density, average stand height, and DBH. Biomass is estimated from a forest stand using stand density multiplied by single-tree biomass. We used allometric equations of dominant species within a forest stand to estimate the single-tree biomass. Field plots co-located with LiDAR footprints were overlaid onto forest stands from the forest inventory data to derive the biomass of these field plots. These field plots were sampled in Tahe in Heilongjiang province (2458 plots within red box in zone I), Lushuihe in Jilin province (6696 plots within red box in zone II), Qinling in Shanxi province (2857 plots within red box in zone III), Xishuangbanna in Yunnan province (393 plots within red box in zone V), and western Sichuan province (10,128 plots within red box in zones IV and VII) from 2003 to 2006.

Table 1. Summary of field data. N_t is the total number of all field forest stands and N_n, N_b, and N_m are the total number of needleleaf forest, broadleaf forest, and mixed forest stands, respectively. Ns is the number of field data listed by source and acquisition time.

Sample Zone	Nt	Nn	Nb	Nm	Biomass Range (Mg·ha−1)	Ns (Source, Acquisition Time)
Cold temperate zone (I)	2415	1506	251	658	8.88–334.10	2355 (NFRI,2003), 60 (field work, 2006)
Temperate zone (II)	6357	1347	3877	1133	7.60–404.08	6307 (NFRI, 2004), 50 (field work, 2007)
Warm temperate zone (III)	2644	332	2005	307	10.60–277.94	2644 (NFRI, 2006)
Subtropical zone (IV)	4475	3609	723	143	9.21–466.67	4475 (NFRI, 2006)
Tropical zone (V)	431	30	401	0	11.16–416.10	391 (NFRI, 2006), 40 (field work, 2010)
Qinghai-Xizang plateau alpine (VII)	6359	3224	2559	576	16.14–509.42	6359 (NFRI, 2006)

Table 1. Summary of field data. N_t is the total number of all field forest stands and N_n, N_b, and N_m are the total number of needleleaf forest, broadleaf forest, and mixed forest stands, respectively. Ns is the number of field data listed by source and acquisition time.

Sample Zone	Nt	Nn	Nb	Nm	Biomass Range (Mg·ha−1)	Ns (Source, Acquisition Time)
Cold temperate zone (I)	2415	1506	251	658	8.88–334.10	2355 (NFRI,2003), 60 (field work, 2006)
Temperate zone (II)	6357	1347	3877	1133	7.60–404.08	6307 (NFRI, 2004), 50 (field work, 2007)
Warm temperate zone (III)	2644	332	2005	307	10.60–277.94	2644 (NFRI, 2006)
Subtropical zone (IV)	4475	3609	723	143	9.21–466.67	4475 (NFRI, 2006)
Tropical zone (V)	431	30	401	0	11.16–416.10	391 (NFRI, 2006), 40 (field work, 2010)
Qinghai-Xizang plateau alpine (VII)	6359	3224	2559	576	16.14–509.42	6359 (NFRI, 2006)

Additional, field measurements were collected in Tahe, Heilongjiang province (zone I) in September 2006; Changbai Mountain, Jilin province (zone II) in June 2007; and Xishuangbanna, Yunnan province (zone V) in October 2010 by the authors’ team. The geographic coordinates of these field measurements were not located in the forest stands from NFRI data to avoid duplicate sampling. Since it was difficult to collect field inventory data across all forested lands in China, field measurements were beneficial as supplementary data. Field measurement schemes that varied by topographical constraints were similar to those used by Lefsky et al. [53] and Nelson et al. [30], but the species, height, and DBH of all trees >5 cm were recorded in each plot. The scheme used in the plots of Tahe and Changbai Mountain employed a 7.5 m radius central plot located on the centroid of a GLAS shot, plus three replicate plots located 20 m away from the center plot at azimuths of 120°, 240°, and 360°. This is appropriate for a plain area with a moderate slope because four plots with a 7.5 m radius in a footprint can be measured easily. Where plots were located on steep slopes in the tropical zone, it was difficult to measure four plots within a footprint, and the scheme employed a single 15 m radius plot located at the center of the GLAS footprint. Although the 15 m radius circle samples only ~21% of the area of a GLAS footprint, the circle’s position is concentrated at the center of the Gaussian distribution of the incident laser beam, and the circle is able to measure trees illuminated by ~76% of the laser. Allometric equations are used to convert field measurements to biomass estimates. In the field plots of Tahe and Changbai Mountain, biomass conversion for most species were based on DBH-only allometric equations [54], except for spruce (Picea koraiensis), which used DBH-H combined equations [55]. Tree heights were difficult to measure because of high-density forest cover and topography constraints in some areas. A sub-plot of spruce was purposefully selected, including high, medium, and low heights to develop the allometric equation (height = 12.423 ln(DBH) − 20.692, R² = 0.818, n = 35) relating height to DBH so that heights could be estimated for all sampled trees, and then average canopy heights of plots could be estimated. Allometric equations used in the plots of Xishuangbanna came from the data sharing system of the Chinese Ecosystem Research Network (http://www.cerndata.ac.cn/). Figure 2 shows the magnitude of aboveground biomass values collected from NFRI data and field measurements across six forest zones. The distribution of AGB values from these data represent an unsystematic and uneven sampling of biomass and cover a large range of AGB values. Biomass values of forest stands in cold temperate zone (Figure 2a) and warm temperate zone (Figure 2c) are normally distributed with closely mean values and standard deviations. In temperate zone (Figure 2b), biomass values are accumulated in the range of (0, 220). The values of biomass in subtropical zone (Figure 2d) are focused in the range from 20 Mg·ha⁻¹ to 120 Mg·ha⁻¹. Forests in Qinghai-Xizang plateau alpine zone (Figure 2f) with AGB in ranges between 20 Mg·ha⁻¹ and 200 Mg·ha⁻¹ have the highest mean biomass and the biggest standard deviation. Biomass values for the tropical zone are comparatively low. There are three aspects that affected the biomass values. First, due to the savage natural environment and logistics constraint in the tropical zone, we were not able to get a large number of field plots or as many as other forest zones. The plots with greatest variability in biomass were difficultly acquired. Second, biomass of rubber plantation is relative low. According to our field measurements, biomass of a 15-year-old stand was less than 60 Mg·ha⁻¹ and ages of most rubber plantation plots were less than 15 years old in the site of field measurements. Third, there is uncertainty in the field survey of NFRI data. There is no dominant species in the region, identification from different investigators could be a potential source of error to estimate biomass.

Figure 2. Distribution of forest aboveground biomass of field inventory data (listed in Table 1) in six forest zones: (a) cold temperate zone, (b) temperate zone, (c) warm temperate zone, (d) subtropical zone, (e) tropical zone, and (f) Qinghai-Xizang plateau alpine zone.

Figure 2. Distribution of forest aboveground biomass of field inventory data (listed in Table 1) in six forest zones: (a) cold temperate zone, (b) temperate zone, (c) warm temperate zone, (d) subtropical zone, (e) tropical zone, and (f) Qinghai-Xizang plateau alpine zone.

2.2.2. LiDAR Waveform Data

To reduce the inconsistency of acquisition time between GLAS observations and MODIS imagery and field inventory data, GLAS data acquired during several periods in 2003 (2A, September to November), 2004 (3A, October to November), 2005 (3C, May to June; 3D, October to November) and 2006 (3F, May to June; 3G, October to November) over China were used in the research. The ICESat program provides 15 data products (GLA01 to GLA15). GLA01 and GLA14 of version 31 were used in the study. GLA01 products provide raw waveforms for each laser shot. The waveforms were recorded in 544 bins with a bin size of 1 ns or 15 cm for the land surface. GLA14 data provide surface elevation, latitude and longitude of footprints, laser range offsets for the signal beginning and end, and location, amplitude, and width of up to six Gaussian peaks. The six Gaussian distributions of a waveform correspond to different vertical structural features of vegetation cover and underlying topography [56,57]. With the record index and shot number, the location (latitude/longitude) of each GLAS footprint from GLA14 data was added to the individual waveform extracted from GLA01 data, as well as other parameters such as noise level and transmitted pulse waveform, which were used later in waveform processing. Usually, winter in northern China lasts approximate six months or longer. Therefore, only GLAS data acquired from May to July (L3C and L3F) were used in forest zones I, II, and VI.

Since waveforms are greatly affected by clouds and system noise, the cloud-contaminated and abnormal waveforms should be identified and discarded before extracting waveform parameters. We estimated the noise levels before the signal beginning and after the signal end from the raw waveform separately using a method based on the histogram (it is a graphical representation of the waveform bins) [56]. Cloudless waveforms were selected using the cloud detection flag (FRir_qaFlag = 15). Saturated signals were identified using the GLAS flag (SatNdx > 0). Then, the waveforms were smoothed using a Gaussian filter with a width similar to the transmitted laser pulse, and the signal beginning and end were identified using a noise threshold [56]. The total waveform energy was calculated by summing all the return energy from the signal beginning to end. Starting from the signal end, the positions of 25% (h25), 50% (h50), and 75% (h75) of accumulated return energy were located [56]. In addition, the heights at which the waveform energy reached 10% (h10), 20% (h20), …, 90% (h90) of total energy were also calculated. Finally, the GLAS data (referenced to TOPEX/Poseidon) were converted to WGS84. In earlier studies, Lefsky et al. [53] chose the maximum canopy height (referred to as h14) squared as the independent variable to develop biomass estimation models. After comparing several GLAS waveform variables, Boudreau et al. [12] concluded that the slope of the leading edge extent (fslope), wflen (defined as waveform length, the offset difference between signal starting and ending [56]), and SRTM (Shuttle Radar Topographic Mission) range were optimal variables for predicting AGB using a regression approach. Sun et al. [42] found a strong relationship between AGB and h50 and h75. These GLAS waveform parameters, along with other GLAS metrics from waveforms that captured the characteristics of the canopy structure used in the study, are given in detail in Table 2.

The vertical extent of each waveform increases as a function of terrain slope and footprint size, along with the change of spatial pattern of ground surface visible to the laser [58]. As the terrain in southern China is generally complex, the waveform parameters usually used in flat areas were not sufficient to obtain biomass estimates with high accuracy in the southern forested regions. Parameters (e.g., leading edge extent [58]; trailing edge extent [58]) influenced by the slope were sensitive to topography. The leading edge extent is defined as the height difference between the elevation of the signal start and the first bin that is at half the maximum signal intensity above noise value, and is related to canopy height variability. The trailing edge extent is determined as the height difference between the signal end and the last bin in which the signal intensity of the waveform is half the maximum intensity. However, when considering waveforms with high differences between canopy and ground signal amplitude, a modified leading and trailing edge extent were appropriate to better represent canopy surface and terrain topography [59]. For example, when regarding waveforms with high ground amplitude and low canopy amplitude, there may be no intersection between the extension of the location at the half maximum amplitude and waveforms of low canopy amplitude. Therefore, the leading edge cannot be calculated. The modified metrics were related to the center peak of the canopy and ground return rather than the location of the half maximum intensity. In addition, a terrain index (h_slope = wflen – 0.5 × 65 × tan(θ)) was applied in this study, defined as the extent difference between waveform length and footprint size (~65 m) multiplied by half the tangent of the slope (θ) (from the 90 m SRTM elevation data) within a footprint [60]. The term terrain index is different from the meaning in Lefsky et al. [53]. The creation of this new variable was based on the idea of equalization treatment for terrain slopes. As mentioned in Sun et al. [56], the maximum canopy height derived from the LiDAR waveform is denoted as the distance from the signal beginning to the ground peak of the waveform (equal to h100). In sloped terrains, the surface slope widens the ground peak of the waveform, and the maximum extent of waveform length can be calculated as the distance from the lowest ground surface to the top of the highest canopy surface in the waveform. This meant that there was a potential offset. Assuming homogeneous canopy heights in the footprint and ground peak at the same height both on the flat surface and sloped terrain, the maximum value of the offset could be calculated as: Offset = 65 × tan(θ). The terrain index is intended to remove the average value of the offset from the waveform length in complex terrain conditions.

2.2.3. MODIS Data

The MODIS (Moderate Resolution Imaging Spectroradiometer) instrument is operating on both the Terra and Aqua satellites. Its detectors measure 36 spectral bands from visible to thermal infrared channels with a viewing swath width of 2230 km and temporal resolution of 1–2 days. One of the MODIS vegetation index products used in this study was MOD13A1, with 16-day image composites at 500 m resolution for 2006, designed to provide consistent, spatial, and temporal comparisons of vegetation conditions that can be used to monitor photosynthetic activity [30,61]. Only the higher quality data, whose nadir-view pixels with minimal residual atmospheric aerosols, cloud-free, were selected for creating composites. The methods of creating composites were introduced by Huete et al. [62]. For pixels with lower quality flags, the best data from the corresponding months over the next two years were selected to fill gaps. The product contains four spectral bands: Blue (459–479 nm), red (620–670 nm), near infrared (841–876 nm) and middle infrared (2105–2155 nm), and two vegetation indices (VIs): Normalized difference vegetation index (NDVI) and enhanced vegetation index (EVI). All images from 23 dates, including early, peak, and leaf-off phenological conditions, and of 19 tiles (h23v04, h23v05, h24v04, h24v05, h25v03, h25v04, h25v05, h25v06, h26v03, h26v04, h26v05, h26v06, h27v04, h27v05, h27v06, h28v05, h28v06, h28v07, and h29v06) were processed to build the entire China land cover satellite dataset. The MODIS Reprojection Tool was used to re-project, mosaic, and resample the data with a pixel size of 500 m.

Two other MODIS data products were also used in this research. One was the vegetation continuous fields (VCF) product (MOD44B) at a spatial resolution of 500 m from 2006. It includes percentage estimates for three land cover types: woody vegetation, herbaceous vegetation, and bare soil [63]. The percentage of wood vegetation cover was used in this study for biomass estimation. Since VCF is a yearly data, the percent of vegetation cover is probably the best expression for continuous distribution of vegetation in space [64,65]. The other MODIS data product was MCD12Q1 for 2006. It provides five global land cover classification systems, which describe land cover properties derived from observations spanning a year’s input of MODIS data at 500 m spatial resolution. As one of the classification systems, the IGBP (International Geosphere Biosphere Programme) global vegetation classification scheme was adopted to organize 17 land cover types in China. All of the land cover types were merged into four classes for the study: needleleaf forests, broadleaf forests, mixed forests and non-forest areas.

2.3. Relating Ground-Based AGB to GLAS Waveform Parameters

All of the GLAS shots in China were filtered as stated above (Section 2.2.2) and masked to remove records that were not in the forest areas. After the selection process, 347,116 GLAS shots remained. GLAS footprints were overlaid onto forest stand maps of biomass derived from field observations to extract observed biomass within the footprints so as to link biomass estimates from ground data to GLAS metrics. Of LiDAR footprints, 80% were randomly selected for training, and the remaining 20% were reserved for verification. Because the effect of forest type on GLAS waveforms, models of biomass estimation from GLAS data were developed separately for needleleaf forests, broadleaf forests and mixed forests in different zones. These models were used to predict biomass for other GLAS footprints that were not co-located with field data.

Stepwise regression was performed to predict AGB in GLAS footprints. Compared with other machine learning methods (such as decision tree or neural networks), the regression method was not limited to the range of data, and showed the ability to yield reasonable results, similar to complex methods [20,66]. To infer GLAS metrics that had high correlation with biomass and to discard parameters that had less effect on biomass estimation, t-tests and F-tests were conducted to identify the predictive waveform parameters and the overall correlation of regression models. Neter and Wasserman [67] described how to use the two methods to make comparisons of regression parameters and test whether any two regression lines are identical. The coefficient of determination (R²) is widely used to evaluation of goodness of fit for the regression models. As the number of explanatory variables increases, the R² values also increase. However, the additional explanatory variables may not be significant, and do not substantially improve R². Therefore, R² is not appropriate for a comparison between models having different numbers of explanatory variables [68]. Adjusted R-square (R²-a) was adopted for the comparisons. The root mean square error (RMSE) was used to assess the predicted AGB versus AGB estimated from field-observations at the footprint level.

To evaluate the performance of biomass estimation from GLAS waveform parameters and assess the impact of terrain on biomass estimation, models of GLAS-derived biomass estimates were developed for forest zones I, II, III, IV, and VII for the following three cases: the first case 1, used GLAS footprints on all slopes and all GLAS waveform parameters (Table 2) except lead, trail, and h_slope; The second case 2, used GLAS footprints on all slopes and all GLAS waveform parameters in Table 2. The goal in comparing the case 1 and case 2 models was to evaluate the behavior of lead, trail, and h_slope, to see whether these metrics were improved biomass estimation in sloped conditions; case 3 split GLAS footprints on SRTM slopes <20° and ≥20° and used all waveform parameters. The purpose of case 3 was to assess the applicability of the split models for biomass estimation by comparing them with models based on case 2. Because most field plots were located on slopes <20° in zone I and on slopes >20° in zone V, GLAS footprints on all slopes were used to train models in these forest zones.

Table 2. Definition of GLAS (the Geoscience Laser Altimeter System) waveform parameters.

GLAS Waveform Parameters	Definition
meanh, medh	Mean and median canopy height calculated from waveform [66].
wflen, centroid, h14	Waveform length, centroid and top tree height from GLA14 product.
fslope, eratio	Front slope and vegetation to surface energy ratio from waveform as calculated by Sun et al. [56]
h10, h20, …, h100	Deciles heights for waveform energy to reach 10%, 20%, …, 100% of total energy starting form signal end [56].
h25, h75	Quartiles heights for waveform energy to reach 25% and 75% of total energy starting form signal end [56].
ht3	Top tree height calculated from waveform with one of corrections as calculated by Sun et al. [56]
gdpamp	Amplitude of Gaussian peaks in GLA14 product.
lead, trail	Leading edge extent and trailing edge extent calculated from waveform [59].
hslope	A terrain index calculated from a function of waveform length, footprints size and slope [60].

Table 2. Definition of GLAS (the Geoscience Laser Altimeter System) waveform parameters.

GLAS Waveform Parameters	Definition
meanh, medh	Mean and median canopy height calculated from waveform [66].
wflen, centroid, h14	Waveform length, centroid and top tree height from GLA14 product.
fslope, eratio	Front slope and vegetation to surface energy ratio from waveform as calculated by Sun et al. [56]
h10, h20, …, h100	Deciles heights for waveform energy to reach 10%, 20%, …, 100% of total energy starting form signal end [56].
h25, h75	Quartiles heights for waveform energy to reach 25% and 75% of total energy starting form signal end [56].
ht3	Top tree height calculated from waveform with one of corrections as calculated by Sun et al. [56]
gdpamp	Amplitude of Gaussian peaks in GLA14 product.
lead, trail	Leading edge extent and trailing edge extent calculated from waveform [59].
hslope	A terrain index calculated from a function of waveform length, footprints size and slope [60].

2.4. Extrapolating GLAS-Derived AGB Estimates to MODIS Imagery

Traditionally, a large number of ground based sample surveys were conducted in a target area and maps of biomass were derived using the forest inventory in combination with land cover maps to assign average value of biomass density to each land cover category. This method may not be adequate for describing the local spatial variability in biomass [48]. Alternatively, a more spatially consistent method to produce biomass maps over large regions was to find the relationships between in situ biomass and a remote sensing signal using a number of statistical models or machine learning techniques and apply the prediction models to each pixel. In our study, we adopted the latter method, to better describe the spatial variability in biomass over China.

The biomass prediction models were applied to all GLAS footprints that were not co-located with field data. AGB estimates from GLAS footprints within a MODIS pixel were averaged as GLAS-derived biomass. Models linking average AGB estimates within MODIS pixels with MODIS imagery reflectance and ancillary data were developed for every forest type in the seven forest zones. Baccini et al. [20] had found that the RMSE decreased as the number of GLAS footprints within a single MODIS pixel increased, and better results for Asian forest were obtained by using MODIS pixels with three GLAS footprints. Therefore, the MODIS-based models were developed using the best case (three GLAS footprints per pixel) in our study. These models were used to predict AGB across all forested MODIS pixels in China. Finally, the forest AGB map was produced at a spatial resolution of 500 m.

Since biomass is closely related to vegetation density, greenness, and structure, vegetation indices associated with biomass were derived from MODIS data. These vegetation indices were developed from both empirical methods such as NDVI and SR (Simple Ratio), and mathematical models e.g., SAVI (Soil Adjusted Vegetation Index). The growing season NDVI, which was determined by the length of the growing season and magnitude of observation, was also an ideal measurement of seasonal greenness [22]. Biomass accumulates during the growing season and reaches its peak at the last stage of the growing season. The growing seasons of zones I, II, and VI were relatively short, from May to August (including seven 16-day periods). In zones III, IV, and VII, growing seasons usually last for 7–8 months (March–October, fifteen 16-day periods). Forests in zone V maintain growing status all year round (equally 23 sensor observations). All the MODIS-derived variables used to develop models for generating the biomass map are listed in Table 3 [69,70,71,72,73,74]. As some of these vegetation indices were functionally redundant in information content [75], a PCA (Principal Components Analysis) approach was adopted to preserve most of the information content of the original images in each 16-day period, and the first three principal components were selected as predictors. This resulted in a total of 21, 45, and 69 PCA transformed predictors for zones I-II-VI, III-IV-VII, and V, respectively. For each MODIS pixel, mean values of reflectance for each band, VCF, SRTM elevation data (aggregated from 90 m to 500 m), and PCA transformed predictors were used to develop Random Forest (RF) regressions. RF has been developed as a new extension of tree-based models in non-parametric statistics and machine learning methods [76,77]. It had been successfully used for biomass estimation in several different contexts [20,76,77,78,79].

Table 3. Vegetation indices derived from MODIS (Moderate Resolution Imaging Spectro-radiometer) data.

Vegetation Indices	Formula	Reference
SR	Band 1/band 2	Colombo et al. [69]
NDVI	(band 2 − band 1)/(band 2 + band 1)	Rouse et al. [70]
EVI	5 × (band 2 – band 1)/(band 2 + 6band 1−7.5band 3 + 1)	Huete et al. [71]
SAVI	1.5 × (band 2 − band 1)/(band 2 + band 1 + 0.5)	Huete. [72]
NDIIb7	(band 2 – band 7)/(band 2 + band 7)	Hunt et al. [73]
ARVI	(band 2 + band 3 − 2band 1)/(band 2 + 2band 1 − band 3)	Kaufman et al. [74]

Table 3. Vegetation indices derived from MODIS (Moderate Resolution Imaging Spectro-radiometer) data.

Vegetation Indices	Formula	Reference
SR	Band 1/band 2	Colombo et al. [69]
NDVI	(band 2 − band 1)/(band 2 + band 1)	Rouse et al. [70]
EVI	5 × (band 2 – band 1)/(band 2 + 6band 1−7.5band 3 + 1)	Huete et al. [71]
SAVI	1.5 × (band 2 − band 1)/(band 2 + band 1 + 0.5)	Huete. [72]
NDIIb7	(band 2 – band 7)/(band 2 + band 7)	Hunt et al. [73]
ARVI	(band 2 + band 3 − 2band 1)/(band 2 + 2band 1 − band 3)	Kaufman et al. [74]

2.5. Accuracy Assessment

The comparison of our national biomass map with forest inventory data from NFRI was performed at two different scales: (1) NFRI data that were not used to train models of GLAS-derived biomass were randomly selected for validation. Because NFRI forest stands had irregular boundaries, forest stands were selected that covered at least two map pixels, to reduce contingency in the comparison. The AGB values of all map pixels within one forest stand were averaged to compare with the observed AGB value of the forest stand. (2) The comparison of national biomass estimates with biomass statistical results was implemented for each of the 31 provinces in China. Biomass estimates of each province were extracted from the Chinese AGB map to compare with biomass statistical data [80] derived from the seventh national forest resource inventory data compiled by the State Forestry Administration of China (2004–2008). The biomass statistical data included below- and aboveground biomass. There is a relationship between aboveground biomass and belowground biomass for a tree of a given species as well as for a given forest or plantation type. The development of belowground biomass to aboveground biomass ratio involves a large human effort and cost. The ratios did not vary significantly with latitude (tropical, temperate or boreal), soil texture (fine, medium or coarse) or tree type (angiosperms or gymnosperms). We used 0.8 to multiply the magnitude of total biomass to calculate AGB according to MacDicken [81]. The flowchart of our study was performed as follows (Figure 3).

Figure 3. Processing flow of national biomass estimation from MODIS, GLAS, and field data.

Figure 3. Processing flow of national biomass estimation from MODIS, GLAS, and field data.

3. Results

3.1. Biomass Estimates from GLAS Waveform Parameters

Logarithm function models were fitted to the needleleaf forests, broadleaf forests and mixed forests in each zone with GLAS waveform parameters as explanatory variables. Table 4 shows the performance of GLAS-based AGB models from case 1 and case 2 (see Section 2.3) in the six forest zones. When GLAS waveform parameters associated with terrain factor (e.g., lead, trail, h_slope) were considered in the regression analysis, the R²-a were significantly increased (numbers in brackets in Table 4). For example, R²-a values increased from 0.503 to 0.591, and from 0.675 to 0.787 in AGB estimation models of needleleaf forests and broadleaf forests, respectively, in the warm temperate zone III. Because of the lack of field survey data and low forest coverage in the Neimeng-Xinjiang arid zone, the models of AGB estimation from GLAS data in this zone were referred to the models in the Qinghai-Xizang plateau alpine zone, which is dominated by spruce and fir trees, similar to the forest types in the Neimeng-Xinjiang arid zone.

To better assess the impact of slope on biomass estimation, the AGB of some forest stands are derived from field survey data in zones II, III, IV, and VII to compare with biomass derived from GLAS data, which are sampled on all slopes, slopes <20° and slopes ≥20° using all GLAS parameters. Average AGB that were estimated from field data in the four forest zones are 134.69 Mg·ha⁻¹, 115.40 Mg·ha⁻¹, 99.16 Mg·ha⁻¹, and 140.24 Mg·ha⁻¹. Slopes are divided into five ranges to illustrate biomass estimate biases on various slopes (Figure 4). In each range, biomass values of all GLAS footprints within one forest stand are compared with the observed biomass value of the forest stand. All forest stands in each slope range are used. AGB estimated from GLAS data on all slopes are close to the biomass estimated from field survey on gentle slopes (0°–20°). When the slope is greater than 20°, the difference between the biomass estimated from GLAS data on all slopes and the field measurement biomass increases. For slopes <20° and ≥20°, the gaps between biomass estimates of model I (low slope), model II (high slope), and field survey biomass are smaller than the gap between biomass estimates of model III (all slopes) and field survey biomass. Moreover, the difference between biomass estimates of model I and field survey biomass is smaller than the difference between biomass estimates of model II and field survey biomass. In general, model III shows a tendency to overestimate AGB with increasing slope. Models I and II show good consistency between AGB estimates and observed biomass in each slope range.

Table 4. R²-a (adjusted R-square) of GLAS-based (the Geoscience Laser Altimeter System) AGB (aboveground biomass) models using datasets of case (1) and case (2) in the six forest zones. The first set of numbers in each group corresponds to the R²-a of AGB models from GLAS footprints on all slopes (without using lead, trail, and h_slope). The second set of numbers (in round brackets) is R²-a values of AGB models from GLAS footprints on all slopes using all GLAS waveform metrics (including lead, trail, and h_slope). NF = needleleaf forest, BF = broadleaf forest, and MF = mixed forest. TMF = tropical monsoon forest. RP = rubber plantation. I: Cold temperate zone, II: Temperate zone, III: Warm temperate zone, IV: Subtropical zone, V: Tropical zone, and VII: Qinghai-Xizang plateau alpine zone.

Forest Zones	R2-a in Case 1 (Case 2)
	NF	BF	MF	RP	TMF
I	0.786 (0.815)	0.658 (0.683)	0.671 (0.689)
II	0.662 (0.688)	0.706 (0.723)	0.694 (0.712)
III	0.503 (0.591)	0.675 (0.787)	0.508 (0.583)
IV	0.552 (0.614)	0.611 (0.660)	0.452 (0.518)
V				0.651 (0.762)	0.624 (0.701)
VII	0.527 (0.599)	0.607 (0.693)	0.557 (0.645)

Table 4. R²-a (adjusted R-square) of GLAS-based (the Geoscience Laser Altimeter System) AGB (aboveground biomass) models using datasets of case (1) and case (2) in the six forest zones. The first set of numbers in each group corresponds to the R²-a of AGB models from GLAS footprints on all slopes (without using lead, trail, and h_slope). The second set of numbers (in round brackets) is R²-a values of AGB models from GLAS footprints on all slopes using all GLAS waveform metrics (including lead, trail, and h_slope). NF = needleleaf forest, BF = broadleaf forest, and MF = mixed forest. TMF = tropical monsoon forest. RP = rubber plantation. I: Cold temperate zone, II: Temperate zone, III: Warm temperate zone, IV: Subtropical zone, V: Tropical zone, and VII: Qinghai-Xizang plateau alpine zone.

Forest Zones	R2-a in Case 1 (Case 2)
	NF	BF	MF	RP	TMF
I	0.786 (0.815)	0.658 (0.683)	0.671 (0.689)
II	0.662 (0.688)	0.706 (0.723)	0.694 (0.712)
III	0.503 (0.591)	0.675 (0.787)	0.508 (0.583)
IV	0.552 (0.614)	0.611 (0.660)	0.452 (0.518)
V				0.651 (0.762)	0.624 (0.701)
VII	0.527 (0.599)	0.607 (0.693)	0.557 (0.645)

AGB estimation models optimized for forests in each of the six forest zones were selected based on the highest adjusted R-square values. The stratified regression equations for AGB estimation from GLAS data are shown and evaluated in Table 5. The results indicate that the adjusted R² of AGB estimation models were greater than 0.68 for the three forest types in most forest zones, exceeding 0.8 for models for needleleaf forests in the cold temperate zone (I), models for broadleaf forests in the warm temperate zone (III) and models for broadleaf forests and mixed forests in the Qinghai-Xizang plateau alpine zone (VII). Values of R²-a range from 0.7 to ~0.8 for all biomass estimation models in the warm temperate zone (II) and tropical zone (V). The adjusted coefficient of determination is lowered to ~0.62 for AGB estimation models of mixed forests in the subtropical zone.

More than 2800 footprints were randomly selected from the validation dataset in the six forest zones. The validation results are shown in Table 6 in terms of forest type. The results show that AGB estimated from GLAS waveforms and biomass converted from forest inventory data were very consistent. R²-a values remained higher than 0.7 for all the forest types in zones I, II, III, and V. RMSE values ranged from 9.38 Mg·ha⁻¹ for needleleaf forests in the cold temperate zone to 43.81 Mg·ha⁻¹ for tropical monsoon rainforests.

Figure 4. Comparison of biomass estimated from models using GLAS footprints on all slopes (purple crosses in model III), slopes <20° (brown squares in model I), and slopes ≥20° (green triangles in model II) with field survey biomass (blue diamonds) in forest zones II (a), III (b), IV (c), and VII (d). Models I and II were trained by GLAS waveform data on slopes <20° and ≥20°, respectively. Models III were developed from GLAS waveform data on all slopes.

Figure 4. Comparison of biomass estimated from models using GLAS footprints on all slopes (purple crosses in model III), slopes <20° (brown squares in model I), and slopes ≥20° (green triangles in model II) with field survey biomass (blue diamonds) in forest zones II (a), III (b), IV (c), and VII (d). Models I and II were trained by GLAS waveform data on slopes <20° and ≥20°, respectively. Models III were developed from GLAS waveform data on all slopes.

Table 5. Optimized models of AGB estimated from GLAS waveform data. R² were the Goodness of Fit for different forest types in the six forest zones. I is cold temperate zone, II is temperate zone, III is warm temperate zone, IV is subtropical zone, V is tropical zone and VII is Qinghai-Xizang plateau alpine zone. (B = biomass, N = number of GLAS footprints sampled by forest plots used in each model).

Zone	Forest Types	Models	R2-a	N
I	Needleaf	ln(B) = 0.7525 medh − 1.6582 meanh + 4.5879 h75 − 2.1283 h80 + 0.4522 h100 + 5.6168	0.821	1203
	Broadleaf	ln(B) = 0.7631 medh + 0.0018 eratio − 0.3965 wflen − 1.8252 h30 + 1.4663 h75 + 2.6246	0.687	237
	Mixed	ln(B) = 0.0096 eratio + 1.2566 h50 − 2.1639 h60 + 1.5212 h70 − 0.4168 h80 + 0.3581 h100 + 3.9187	0.702	548
II	Needleaf	ln(B) = 0.6035 meanh + 0.0856 trail – 0.9138 h20 + 1.0116 h50 + 0.2132 h100 + 4.8012 (slope <20°)	0.711	716
		ln(B) = 0.7105 medh + 0.1521 lead – 0.1292 trail – 0.7469 h25 + 1.2734 h50 + 0.2516 h100 + 4.6931 (slope ≥20°)	0.718	542
	Broadleaf	ln(B) = 0.3427 lead + 0.1528 h100 − 0.0637 eratio + 0.3829 h50 − 0.6614 h10 − 1.2331 ln(h100) + 4.5864 (slope <20°)	0.784	1920
		ln(B) = 0.3861 lead + 0.3317 trail + 0.1592 h100 – 0.1173 eratio + 0.1243 h50 – 0.3724 h20 – 0.8618 ln(h100) + 4.2077 (slope ≥20°)	0.755	1566
	Mixed	ln(B) = 0.0842 h75 − 0.1732 wflen − 0.1574 h25 + 0.2066 h100 + 4.7751 (slope <20°)	0.756	624
		ln(B) = 0.0614 h75 + 0.3059 lead – 0.1028 trail – 0.0952 h25 + 0.1337 h100 + 4.6832 (slope ≥20°)	0.737	226
II	Needleaf	ln(B) = 0.0683 medh − 0.0355 lead − 0.0168 trail − 3.3667 wflen + 7.3146 h100 − 3.3629 ht3 − 0.0071 hslope + 0.4835 (slope <20°)	0.694	133
		ln(B) = 0.0588 meanh – 0.0517 lead – 0.0269 trail + 8.5293 h100 – 2.8243 wflen – 5.9678 h50 – 0.0124 hslope + 0.5139 (slope ≥20°)	0.701	125
	Broadleaf	ln(B) = 0.0082 meanh − 0.0058 medh + 0.0019 centroid + 0.0008 eratio − 0.0039 h50 + 0.0359 h75 + 0.0021 hslope − 0.0005 (ht3)2 + 2.9672 (slope <20°)	0.898	1116
		ln(B) = 0.0078 medh + 0.0645 lead − 0.0616 trail − 0.0069 h50 + 0.0022 eratio + 0.0295 h75 + 0.0016 hslope – 0.0004 (ht3)2 + 2.658 (slope ≥20°)	0.884	696
	Mixed	ln(B) =0.2178 wflen + 0.1335 h50 − 0.2532 h14 + 0.0084 gdpamp − 0.0010 (h75)2 + 0.0006 (h100)2 − 0.0006 hslope + 4.7153 (slope <20°)	0.703	128
		ln(B) = 0.1832 wflen + 0.1049 h50 − 0.2557 h14 + 0.0078 gdpamp − 0.0016 (h75)2 + 0.0005 (h100)2 − 0.0008 hslope + 4.3551 (slope ≥20°)	0.681	117
IV	Needleaf	ln(B) = 0.0778 h14 − 0.0706 trail – 0.0436 h50 + 0.0551 h100 − 0.0642 hslope + 0.0010 (h75)2 + 2.3814 (slope < 20°)	0.716	1582
		ln(B) = 0.0641wflen + 0.0137lead − 0.0626trail + 0.0396h100 − 0.0845hslope + 0.0013 (h75)2 + 2.4783 (slope ≥20°)	0.707	1133
	Broadleaf	ln(B) = 0.0375 wflen + 0.0584 centroid − 0.0263 fslope + 0.0625 h100 − 0.0563 ht3 − 0.0405 hslope + 3.9661 (slope <20°)	0.715	258
		ln(B) = 0.0406 h14 + 0.0818 lead − 0.0282 trail − 0.0249 fslope + 0.0588 h100 − 0.0571 ht3 − 0.0331 hslope + 3.9512 (slope ≥20°)	0.703	212
	Mixed	ln(B) = 0.0821 lead + 0.0563 fslope − 0.0903 trail + 0.0426 h100 − 0.0022 hslope + 4.3916 (slope <20°)	0.624	215
		ln(B) = 0.0865 lead + 0.0556 fslope − 0.0941 trail + 0.0483 h100 − 0.0029 hslope + 4.7523 (slope ≥20°)	0.646	184
V	Tropical monsoon forest	ln(B) = 0.1856 meanh + 0.0166 eratio + 0.1043 h100 − 0.096 h14 − 0.0279 lead − 0.1345 h25 + 2.3628	0.759	202
	Rubber plantation	ln(B) = 0.8636h50 – 0.7568trail + 2.8519	0.708	75
VII	Needleaf	ln(B) = 0.0525 centroid − 0.0232 wflen + 0.0221 hslope + 3.2316 (slope <20°)	0.751	1303
		ln(B) = 0.0925 lead − 0.0252 trail − 0.0237 h14 + 0.0242 hslope + 3.3106 (slope ≥20°)	0.749	975
	Broadleaf	ln(B) = 2.1091 meanh − 0.5753 medh + 0.1024 centroid − 1.4145 h50 + 0.3168 h100 − 0.071 ht3 − 0.040 hslope + 2.3511 (slope <20°)	0.852	1011
		ln(B) = 2.0882 meanh − 0.5249 medh − 1.3872 h50 + 0.3455 h100 − 0.036 ht3 − 0.037 hslope + 2.2859 (slope ≥20°)	0.846	948
	Mixed	ln(B) = 0.0527 meanh + 0.1237 medh + 2.6855 (slope <20°)	0.811	188
		ln(B) = 0.0487 meanh + 0.1025 medh + 0.0536 lead – 0.0174 trail + 2.6685 (slope ≥20°)	0.802	181

Table 5. Optimized models of AGB estimated from GLAS waveform data. R² were the Goodness of Fit for different forest types in the six forest zones. I is cold temperate zone, II is temperate zone, III is warm temperate zone, IV is subtropical zone, V is tropical zone and VII is Qinghai-Xizang plateau alpine zone. (B = biomass, N = number of GLAS footprints sampled by forest plots used in each model).

Zone	Forest Types	Models	R2-a	N
I	Needleaf	ln(B) = 0.7525 medh − 1.6582 meanh + 4.5879 h75 − 2.1283 h80 + 0.4522 h100 + 5.6168	0.821	1203
	Broadleaf	ln(B) = 0.7631 medh + 0.0018 eratio − 0.3965 wflen − 1.8252 h30 + 1.4663 h75 + 2.6246	0.687	237
	Mixed	ln(B) = 0.0096 eratio + 1.2566 h50 − 2.1639 h60 + 1.5212 h70 − 0.4168 h80 + 0.3581 h100 + 3.9187	0.702	548
II	Needleaf	ln(B) = 0.6035 meanh + 0.0856 trail – 0.9138 h20 + 1.0116 h50 + 0.2132 h100 + 4.8012 (slope <20°)	0.711	716
		ln(B) = 0.7105 medh + 0.1521 lead – 0.1292 trail – 0.7469 h25 + 1.2734 h50 + 0.2516 h100 + 4.6931 (slope ≥20°)	0.718	542
	Broadleaf	ln(B) = 0.3427 lead + 0.1528 h100 − 0.0637 eratio + 0.3829 h50 − 0.6614 h10 − 1.2331 ln(h100) + 4.5864 (slope <20°)	0.784	1920
		ln(B) = 0.3861 lead + 0.3317 trail + 0.1592 h100 – 0.1173 eratio + 0.1243 h50 – 0.3724 h20 – 0.8618 ln(h100) + 4.2077 (slope ≥20°)	0.755	1566
	Mixed	ln(B) = 0.0842 h75 − 0.1732 wflen − 0.1574 h25 + 0.2066 h100 + 4.7751 (slope <20°)	0.756	624
		ln(B) = 0.0614 h75 + 0.3059 lead – 0.1028 trail – 0.0952 h25 + 0.1337 h100 + 4.6832 (slope ≥20°)	0.737	226
II	Needleaf	ln(B) = 0.0683 medh − 0.0355 lead − 0.0168 trail − 3.3667 wflen + 7.3146 h100 − 3.3629 ht3 − 0.0071 hslope + 0.4835 (slope <20°)	0.694	133
		ln(B) = 0.0588 meanh – 0.0517 lead – 0.0269 trail + 8.5293 h100 – 2.8243 wflen – 5.9678 h50 – 0.0124 hslope + 0.5139 (slope ≥20°)	0.701	125
	Broadleaf	ln(B) = 0.0082 meanh − 0.0058 medh + 0.0019 centroid + 0.0008 eratio − 0.0039 h50 + 0.0359 h75 + 0.0021 hslope − 0.0005 (ht3)2 + 2.9672 (slope <20°)	0.898	1116
		ln(B) = 0.0078 medh + 0.0645 lead − 0.0616 trail − 0.0069 h50 + 0.0022 eratio + 0.0295 h75 + 0.0016 hslope – 0.0004 (ht3)2 + 2.658 (slope ≥20°)	0.884	696
	Mixed	ln(B) =0.2178 wflen + 0.1335 h50 − 0.2532 h14 + 0.0084 gdpamp − 0.0010 (h75)2 + 0.0006 (h100)2 − 0.0006 hslope + 4.7153 (slope <20°)	0.703	128
		ln(B) = 0.1832 wflen + 0.1049 h50 − 0.2557 h14 + 0.0078 gdpamp − 0.0016 (h75)2 + 0.0005 (h100)2 − 0.0008 hslope + 4.3551 (slope ≥20°)	0.681	117
IV	Needleaf	ln(B) = 0.0778 h14 − 0.0706 trail – 0.0436 h50 + 0.0551 h100 − 0.0642 hslope + 0.0010 (h75)2 + 2.3814 (slope < 20°)	0.716	1582
		ln(B) = 0.0641wflen + 0.0137lead − 0.0626trail + 0.0396h100 − 0.0845hslope + 0.0013 (h75)2 + 2.4783 (slope ≥20°)	0.707	1133
	Broadleaf	ln(B) = 0.0375 wflen + 0.0584 centroid − 0.0263 fslope + 0.0625 h100 − 0.0563 ht3 − 0.0405 hslope + 3.9661 (slope <20°)	0.715	258
		ln(B) = 0.0406 h14 + 0.0818 lead − 0.0282 trail − 0.0249 fslope + 0.0588 h100 − 0.0571 ht3 − 0.0331 hslope + 3.9512 (slope ≥20°)	0.703	212
	Mixed	ln(B) = 0.0821 lead + 0.0563 fslope − 0.0903 trail + 0.0426 h100 − 0.0022 hslope + 4.3916 (slope <20°)	0.624	215
		ln(B) = 0.0865 lead + 0.0556 fslope − 0.0941 trail + 0.0483 h100 − 0.0029 hslope + 4.7523 (slope ≥20°)	0.646	184
V	Tropical monsoon forest	ln(B) = 0.1856 meanh + 0.0166 eratio + 0.1043 h100 − 0.096 h14 − 0.0279 lead − 0.1345 h25 + 2.3628	0.759	202
	Rubber plantation	ln(B) = 0.8636h50 – 0.7568trail + 2.8519	0.708	75
VII	Needleaf	ln(B) = 0.0525 centroid − 0.0232 wflen + 0.0221 hslope + 3.2316 (slope <20°)	0.751	1303
		ln(B) = 0.0925 lead − 0.0252 trail − 0.0237 h14 + 0.0242 hslope + 3.3106 (slope ≥20°)	0.749	975
	Broadleaf	ln(B) = 2.1091 meanh − 0.5753 medh + 0.1024 centroid − 1.4145 h50 + 0.3168 h100 − 0.071 ht3 − 0.040 hslope + 2.3511 (slope <20°)	0.852	1011
		ln(B) = 2.0882 meanh − 0.5249 medh − 1.3872 h50 + 0.3455 h100 − 0.036 ht3 − 0.037 hslope + 2.2859 (slope ≥20°)	0.846	948
	Mixed	ln(B) = 0.0527 meanh + 0.1237 medh + 2.6855 (slope <20°)	0.811	188
		ln(B) = 0.0487 meanh + 0.1025 medh + 0.0536 lead – 0.0174 trail + 2.6685 (slope ≥20°)	0.802	181

Table 6. AGB predicted from GLAS data were compared with biomass converted from forest inventory data in six forest zones. R² and RMSE (in parentheses) values correspond to each forest type within six forest zones. NF = needleleaf forest, BF = broadleaf forest, and MF = mixed forest, RP = rubber plantation, TMF = tropical monsoon forest. Numbers (R² and RMSE) with superscripts 1 and 2 correspond to estimation models from GLAS footprints on slopes <20° and ≥20°. I is cold temperate zone, II is temperate zone, III is warm temperate zone, IV is subtropical zone, V is tropical zone and VII is Qinghai-Xizang plateau alpine zone.

Forest Zones	R2 and RMSE (Mg·ha−1) in AGB Models Validation
	NF	BF	MF	RP	TMF
I	0.861 (9.38)	0.854 (12.92)	0.864 (13)
II	0.798 (15.64) 1	0.837 (12.94) 1	0.901 (11.25) 1
	0.786 (20.82) 2	0.806 (19.79) 2	0.869 (18.18) 2
III	0.791 (12.47) 1	0.725 (14.26) 1	0.808 (9.28) 1
	0.783 (18.56) 2	0.739 (20.64) 2	0.798(22.87) 2
IV	0.725 (23.60) 1	0.706 (21.91) 1	0.642 (25.57) 1
	0.737 (21.85) 2	0.696 (20.18) 2	0.635 (23.76)2
V				0.887 (12.08)	0.807 (43.81)
VII	0.712 (18.76) 1	0.679 (25.63) 1	0.705 (24.08) 1
	0.718 (20.63)	0.682 (26.04)	0.682 (25.27)

Table 6. AGB predicted from GLAS data were compared with biomass converted from forest inventory data in six forest zones. R² and RMSE (in parentheses) values correspond to each forest type within six forest zones. NF = needleleaf forest, BF = broadleaf forest, and MF = mixed forest, RP = rubber plantation, TMF = tropical monsoon forest. Numbers (R² and RMSE) with superscripts 1 and 2 correspond to estimation models from GLAS footprints on slopes <20° and ≥20°. I is cold temperate zone, II is temperate zone, III is warm temperate zone, IV is subtropical zone, V is tropical zone and VII is Qinghai-Xizang plateau alpine zone.

Forest Zones	R2 and RMSE (Mg·ha−1) in AGB Models Validation
	NF	BF	MF	RP	TMF
I	0.861 (9.38)	0.854 (12.92)	0.864 (13)
II	0.798 (15.64) 1	0.837 (12.94) 1	0.901 (11.25) 1
	0.786 (20.82) 2	0.806 (19.79) 2	0.869 (18.18) 2
III	0.791 (12.47) 1	0.725 (14.26) 1	0.808 (9.28) 1
	0.783 (18.56) 2	0.739 (20.64) 2	0.798(22.87) 2
IV	0.725 (23.60) 1	0.706 (21.91) 1	0.642 (25.57) 1
	0.737 (21.85) 2	0.696 (20.18) 2	0.635 (23.76)2
V				0.887 (12.08)	0.807 (43.81)
VII	0.712 (18.76) 1	0.679 (25.63) 1	0.705 (24.08) 1
	0.718 (20.63)	0.682 (26.04)	0.682 (25.27)

3.2. National Forest AGB Estimation

GLAS-derived estimates of biomass were used to train regression models that were used to generate a national biomass map from MODIS imagery and VCF using the RF algorithm. To better show the general spatial characteristics of national AGB distribution, the biomass map is converted into eight ranges (Figure 5). The four provinces with the highest AGB are Yunnan, Heilongjiang, Xizang, and Sichuan. The total AGB of each of these four provinces exceeds 1000 Mt (million tons), accounting for 13.2%, 11.7%, 10.2%, and 9.7% of the national total AGB, respectively. With respect to average AGB (Mg·ha⁻¹), Xizang, Hainan, and Jilin are the top three provinces, with values of 139.26 Mg·ha⁻¹, 121.53 Mg·ha⁻¹, and 114.91 Mg·ha⁻¹, respectively.

The AGB of most forests in China ranges between 60.1 Mg·ha⁻¹ and 90 Mg·ha⁻¹ (Figure 6a), accounting for 37.8% of the total national AGB. In forest zones I, III, and IV, biomass in this range account for 56.8%, 61.2%, and 49.5% of total biomass, respectively (Figure 6b). The second largest biomass is in the range from 90.1 Mg·ha⁻¹ to 120 Mg·ha⁻¹, nearly a quarter (23.6%) of the total Chinese AGB. The magnitude of biomass in this range represents the largest proportion of total biomass in forest zone II. Forests with AGB greater than 90 Mg·ha⁻¹ hold more than half of the total national biomass. Forest AGB ranging from 120.1 Mg·ha⁻¹ to 150 Mg·ha⁻¹ and from 150.1·Mg·ha⁻¹ to 200 Mg·ha⁻¹ are at the same magnitude, in that the total AGB values for these two ranges were between 1600 Mt and 1850 Mt. The same applies for AGB ranges of 30.1 Mg·ha⁻¹ to 60 Mg·ha⁻¹ and 200.1 Mg·ha⁻¹ to 300 Mg·ha⁻¹, with the total AGB value for both ranging between 750 Mt and 850 Mt. The total biomass in the highest (biomass greater than 300 Mg·ha⁻¹) and lowest (biomass less than 30 Mg·ha⁻¹) AGB classes account for only a small part of the total AGB. Forests with AGB greater than 200 Mg·ha⁻¹ include 5.9% of the total Chinese biomass, and cover some areas of zones I and II, southwestern Sichuan province in zone IV, southern Yunnan province in zone V, and southeast of Xizang in zone VII. In terms of biomass greater than 150 Mg·ha⁻¹, zones V, VII, and IV are the top three zones.

Figure 5. Map of aboveground biomass of Chinese forest. (The word “sheng” is spelled in Chinese Pinyin, and means “province” in English; similar to “shi” and “zizhiqu”).

Figure 5. Map of aboveground biomass of Chinese forest. (The word “sheng” is spelled in Chinese Pinyin, and means “province” in English; similar to “shi” and “zizhiqu”).

Figure 6. Distribution of forest AGB in seven forest zones in eight classes.

Figure 6. Distribution of forest AGB in seven forest zones in eight classes.

3.3. Assessment of National AGB

More than 3000 forest stands derived from NFRI data in four forest zones (zone I, II, III, and VII) with AGB values ranging from ~30 Mg·ha⁻¹ to ~270 Mg·ha⁻¹ are used to validate the biomass map derived from GLAS-derived biomass and MODIS data (Figure 7). As it is possible that more than one map pixel falls into one forest stand, the biomass value of a forest stand was compared with the averaged AGB value of two or more pixels. High correlation was found for Lushuihe, Jilin province, with an R² of 0.672, and the minimum RMSE is seen in Tahe, Heilongjiang province, with a value of 8.05 Mg·ha⁻¹. The total AGB summed for each province (including municipalities directly under the central government) are compared with statistical biomass derived from forest inventory data. R-square with 0.9429 shows overall high consistency (Figure 8). The estimate of total AGB for Chinese forests around 2006 was about 12,622 Mt, almost equal to the statistical results (12,617 Mt) of the seventh national forest inventory data (2004–2008). Fifteen provinces had relative errors (percentage) less than 20%, and the relative errors of ten provinces were between 20% and 50%, indicating overall high consistency between estimated AGB and statistical results at the province scale (Table 7). The other six provinces had relative errors greater than 50%, but less than 100%.

Figure 7. Comparison of forest inventory biomass (measured biomass) with GLAS- and MODIS-derived national map of AGB (predicted biomass) in four forest zones. (a) Tahe in zone I; (b) Luishuihe in zone II; (c) Qinling in zone III; (d) Western Sichuan in zone VII.

Figure 7. Comparison of forest inventory biomass (measured biomass) with GLAS- and MODIS-derived national map of AGB (predicted biomass) in four forest zones. (a) Tahe in zone I; (b) Luishuihe in zone II; (c) Qinling in zone III; (d) Western Sichuan in zone VII.

Figure 8. Comparison of biomass derived from field data with AGB prediction from MODIS and GLAS data using biomass of 31 provinces.

Figure 8. Comparison of biomass derived from field data with AGB prediction from MODIS and GLAS data using biomass of 31 provinces.

Table 7. Comparison of estimated total AGB form national map with statistical biomass in 31 provinces. Relative error is the ratio of AGB difference of estimates and forest inventory data to the biomass of forest inventory data. (Due to lack of forest inventory data in Taiwan, Hong Kong and Macau, AGB of these regions were not listed).

Province	Relative Error (%)	Province	Relative Error (%)	Province	Relative Error (%)
Anhui	6.38	Heilongjiang	−2.86	Qinghai	−68.53
Beijing	8.08	Henan	−48.29	Shaanxi	13.56
Chongqing	11.13	Hubei	76.63	Shandong	−88.21
Fujian	6.24	Hunan	22.07	Shanghai	−82.08
Gansu	−0.77	Jiangsu	−70.18	Shanxi	−13.03
Guangdong	19.65	Jiangxi	7.43	Sichuan	5.72
Guangxi	3.89	Jilin	3.75	Tianjin	−23.72
Guizhou	−15.12	Liaoning	49.21	Xinjiang	−31.07
Hainan	32.56	Neimenggu	−28.35	Xizang	−15.45
Hebei	−27.73	Ningxia	−48.37	Yunnan	21.95
				Zhejiang	51.84

Table 7. Comparison of estimated total AGB form national map with statistical biomass in 31 provinces. Relative error is the ratio of AGB difference of estimates and forest inventory data to the biomass of forest inventory data. (Due to lack of forest inventory data in Taiwan, Hong Kong and Macau, AGB of these regions were not listed).

Province	Relative Error (%)	Province	Relative Error (%)	Province	Relative Error (%)
Anhui	6.38	Heilongjiang	−2.86	Qinghai	−68.53
Beijing	8.08	Henan	−48.29	Shaanxi	13.56
Chongqing	11.13	Hubei	76.63	Shandong	−88.21
Fujian	6.24	Hunan	22.07	Shanghai	−82.08
Gansu	−0.77	Jiangsu	−70.18	Shanxi	−13.03
Guangdong	19.65	Jiangxi	7.43	Sichuan	5.72
Guangxi	3.89	Jilin	3.75	Tianjin	−23.72
Guizhou	−15.12	Liaoning	49.21	Xinjiang	−31.07
Hainan	32.56	Neimenggu	−28.35	Xizang	−15.45
Hebei	−27.73	Ningxia	−48.37	Yunnan	21.95
				Zhejiang	51.84

4. Discussion

A method to assess the AGB of Chinese forests was developed by combining large footprint LiDAR data, field measurements, and optical remote sensing data. The results demonstrated that using a small part of forest inventory data to estimate national AGB was feasible, and the estimates were reasonable. GLAS footprint data as a method of sampling were important for linking field plot data to MODIS data because the acquisition of field plot data is time-consuming, and LiDAR samples could extrapolate biomass of limited field plots to more footprints in a greater range. The LiDAR waveform records the interaction of a laser beam with trees in a footprint, starting from the crown. The crown shape reflects the basic geometric structure of a tree. The crown shape of a needleleaf tree is similar to a cone. Therefore, to estimate biomass within a GLAS footprint located in a needleleaf forest, fewer waveform parameters reflecting the characteristics of canopy vertical structure were used in the regression model, such as the mean canopy height or h50. The structure of a broadleaf tree is more complex, and cannot be described by one type of tree crown. More parameters reflecting the overall structure of canopies were used in the models. Moreover, we found waveform energy heights of 75% (h75) were positively correlated with biomass, and used them in the estimation models. Other researchers have also found that h75 was a good indicator of forest vertical structure [46]. In addition, a new waveform parameter h_slope proved to be helpful for predicting footprint-level biomass. In study of boreal forest biomass with large footprint LiDAR at similar latitudes, Boudreau et al. [12] used variables of GLAS waveform extent, waveform front rising slope angle, and a terrain index to estimate biomass in Quebec, Canada. We determined the mean height, median height, and h100 (similar to waveform extent) as the main parameters for predicting AGB in zone I at the same latitude as Quebec. Lefsky et al. [36,53] found that the mean canopy height and the square of the median canopy height had high correlations with AGB in Oregon, USA. We discovered that h50 (similar to median height) and h75 were useful for estimating AGB in forest zone II, at a latitude similar that of Oregon.

The GLAS models tend to overestimate field aboveground biomass with increasing slope [56,58]. The GLAS-sampled biomass was divided into three classes according to the slope in the footprints (all slopes, slopes <20° and slopes ≥20°) to analyze the impact of slopes on AGB estimates. Figure 4 shows that AGB estimates using footprints across all slopes compared well with biomass derived from field data on moderate slopes (0°–20°). This fact may partly be due to logging activity in areas with low to moderate slopes, leading to younger, secondary forests with coverage condition simplifying estimation of biomass. AGB estimated from GLAS footprints on slopes <20° was closer to the biomass derived from field data than all other slopes (Figure 4). In steeply sloped areas, better results were achieved using GLAS footprints on slopes ≥20° rather than on all slopes. This could mean that slopes on <20° and ≥20° combined are more appropriate for predicting biomass than using footprints across all slopes. Another issue that should be noted is that, the biomass range of field data may not to include high values of AGB under moderate terrain conditions of slope <20°, and could be a potential source of error in estimating AGB of GLAS footprints. This uncertainty should be checked before training the AGB estimation models to ensure the representativeness of the training samples.

After calculating biomass from the GLAS footprints, the GLAS-based estimates of biomass can serve as a calibrated inventory dataset for all of China, overcoming the dependency on national forest inventory data. Extrapolating a large number of GLAS-based AGB estimates to a continuous map spanning China requires widespread surface coverage. MODIS image products with coarser spatial resolution and high spectral resolution were useful for scaling forest AGB from plots to sub-continental scales. MODIS images used in this study included optical bands centered on visible, near infrared, and middle infrared, that many studies have indicated are sensitive to forest cover and forest stand structure. MODIS images in high forest cover areas saturate at relatively low levels of forest biomass. Fortunately, good results in AGB estimation of GLAS footprints could offer reliable biomass samples for training across the full range of AGB within a MODIS image. This could improve the saturation of reflectance values in densely forested areas more than using only MODIS images and forest inventory data for biomass estimation.

As expected, the highest forest AGB was found in the southwest of Xizang province (>300 Mg·ha⁻¹), while the lowest forest AGB was found in Hebei province in the north, and Jiangsu province in the south (<50 Mg·ha⁻¹). Good agreement (relative error <20%) between estimates of AGB and forest inventory data was seen in most southern provinces of China. Large discrepancies mainly occurred in provinces located in north-central and east-central China. The potential reason for this is low forest coverage in these areas. The low cover could lead to mistakes in forest/non-forest classification, producing significant differences in AGB estimation. As shown in Figure 6, the biomass ranges were different in seven forest zones. In zone I, the Northern Da Hinggan Mountains suffered from serious fires in 1987. Nearly 20 years later, a field survey was conducted in 2006, and showed that regenerated forest coverage had improved significantly under intensive management, but the stand density and canopy height had not reached levels that existed before the big fire. Before the big fire, the average biomass of larch (Larix gmelinii (Rupr.) Kuzen.) in this area was about 105.05 Mg·ha⁻¹ [55]. After the big fire, it is expected that forest biomass density would decreased. In our research, a biomass range between 60.1 Mg·ha⁻¹ and 90 Mg·ha⁻¹ accounted for the major part of total AGB, consistent with the result of a previous study that the average biomass of a natural young larch forest was about 66.93 Mg·ha⁻¹ [55]. The Xiao Hinggan Mountains has not undergone any large natural disasters in the last 20 years. Thus, forests in this region relatively undestroyed compared with the Da Hinggan Mountains. Therefore, the largest biomass range, from 90.1 Mg·ha⁻¹ to 120 Mg·ha⁻¹, occurs in zone II. The second largest biomass range was from 120.1 to 150 Mg·ha⁻¹. The composition of biomass ranges in zone III was similar to that in zone I. Because the area of forested land in zone III was smaller than in zone I, the amount of total AGB in zone III was less than in zone I. Affected by the long-term human activity, forest zone III was dominated by secondary and man-made forests. Chinese pine (Pinus tabulaeformis Carr.), armand pine (Pinus armandi Franch.), and Japanese red pine (Pinus densiflora Sieb. et Zucc.) were the major species of needleleaf trees, while typical species of broadleaf trees were oak (Quercus L.) and birch (Betula L.). The average biomass of these species at the age of about 30 years were 96.37 Mg·ha⁻¹, 88.9 Mg·ha⁻¹, 89.27 Mg·ha⁻¹ (at the age of 60 years), 49.13 Mg·ha⁻¹, and 84.47 Mg·ha⁻¹, respectively [55]. This is consistent with our conclusion that most of the forest biomass was in the range between 60 Mg·ha⁻¹ and 90 Mg·ha⁻¹. Forest zone IV had the largest forested area of all forest zones, and the largest total amount of biomass, ranging from 60.1 Mg·ha⁻¹ to 90 Mg·ha⁻¹. Of the seven forest zones, zone IV had the highest total biomass because it was dominated by middle-aged man-made forests and young secondary forests. Biomass ranges between 150.1 Mg·ha⁻¹ to 200 Mg·ha⁻¹ accounted for the greatest amount of biomass in zones V and VII. Zone V was mainly covered by tropical monsoon forests with high stand density and canopy height. The forest region in zone VII was located in the southeastern Xizang province, and retains large areas of virgin forest. The altitude of the mountain valley is relative low, and the region has high levels of precipitation and accumulated temperature. All of these conditions result in a high biomass range (>150 Mg·ha⁻¹) in zone VII.

It was a difficult to evaluate the national forest AGB. In order to comprehensively evaluate the AGB estimates, validation was carried out at two scales, the forest stand and nationwide scales. As much as possible, forest stands with the greatest ranges in spatial distribution of biomass were selected. All values of R² were lower than 0.7, mainly because more than one pixel were dropped in one forest stand, so the value of biomass in the stand was compared with the average AGB value of two or more pixels. The R² values were below 0.55 in the Tahe and Qinling regions. The low R² value in Tahe was mainly due to mature forest suffered great loss during the big fire in 1987. Relatively low stand density and canopy heights would affect biomass estimates in GLAS footprints. The field plots in the Qinling area were located on the north slope of Qinling Mountain. The complex topography makes it difficult to estimate AGB with high accuracy. However, the variances of the forest stands in Tahe and Qinling were lower than for that in Lushuihe and western Sichuan. It may be inferred that forest stands in Tahe and Qinling were more homogeneous than that of the other two places. There were small errors between the national AGB estimation and forest inventory data (12,622 Mt vs. 12,617 Mt). A total of 415,000 million fixed forest inventory sample plots (http://www.forestry.gov.cn/) were measured in the seventh forest resource investigation, while fewer than 30,000 plots were used in this research, with the intent of showing that a reasonable and reliable AGB map can be predicted by a limited number of field plots and satellite data. Indeed, large differences existed in some provinces, and more studies are needed at the next step, with more detailed forest inventory data for each province. The forest AGB map would be useful for depicting and quantifying the distribution of forest aboveground biomass over entire landscapes in China.

To reduce errors, LiDAR data were selected during the growing seasons, according to geographical forest distribution. For example, in high-latitude or cold areas, a summer GLAS dataset was selected. In contrast, at mid-latitudes, all datasets except for winter datasets were used for modeling. In the process of extrapolating from GLAS biomass samples to the 500-m spatial scale, MODIS images were processed in a similar manner to maintain the consistency of the acquisition phase. There are some other limitations in the forest AGB estimation. First, Chinese forested lands were roughly divided into seven zones. There are more than 15 eco-regions associated with forests in the Chinese Vegetation Atlas [49]. Field measurements in these regions were time-consuming and labor-intensive. Regions aggregated in several larger zones were required for macroscopic study using remote sensing techniques. We aggregated the forests in China into seven zones based on the geographic characteristics of the forested areas [49,51] to ensure the rationality of partitioning. Second, the forest inventory data represented an uneven sampling of spatial distribution because of human resource constraints and data access restrictions. Third, the acquisition time of GLAS data, MODIS images, and field plots were different in 2005 and 2006, and from 2003 to 2010 (most were between 2003 and 2007, except the field survey in Yunnan in 2010). However, most field data were collected at the end of the growing season between 2003 and 2007 to obtain the biomass at the peak of accumulation, so this could not be a major source of the error. Therefore, seasonal differences in acquisition time between GLAS data and MODIS imagery was likely a potential bias. As GLAS data were acquired over a broad range of latitudes, longitudes, and elevations, the effect of the presence or absence of leaves in GLAS data had to be considered for generating AGB map. Furthermore, snow might be present in high-latitude or high-elevation areas, which could introduce errors in the regression between GLAS-derived biomass and MODIS data.

5. Conclusions

We used GLAS waveform data and MODIS imagery to generate an accurate forest AGB map of China. A statistical relationship between field biomass estimates and GLAS metrics was examined. When GLAS waveform parameters associated with terrain factors (e.g., lead, trail, h_slope) were considered in a regression analysis, the R²-a of the estimation models were significantly increased (e.g., R²-a changed from 0.698 to 0.805 in biomass estimation model of broadleaf forests in zone III). In biomass estimation from GLAS waveform data on slopes <20° and ≥20°, the bias between biomass estimates of the models and biomass estimated from field data were smaller than the bias between biomass estimates of models using GLAS footprints on all slopes and biomass estimated from field data. The assessment of the national biomass map using forest inventory data was performed at two different aggregated scales: (1) At the forest stand level, high correlation appeared in Lushuihe, Jilin province with an R² of 0.672, and the smallest RMSE was found in Tahe, Heilongjiang province, with a value of 8.05 Mg·ha⁻¹. (2) Comparison of national biomass estimates summed in each province with statistical biomass derived from forest inventory data was performed for 31 provinces. Fifteen provinces had relative errors less than 20%, and relative errors of ten provinces were between 20% and 50%. The other six provinces had relative errors greater than 50%, but less than 100%. These result showed overall high consistency between estimated AGB and statistical results at the province scale.

In summary, this study demonstrated that GLAS and MODIS data can be used to estimate forest aboveground biomass at a national scale. GLAS footprints as a way of sampling can be extrapolated to estimate biomass for all of China at a spatial resolution of 500 m. The plot-based methodology using limited ground measurements would reduce the dependence on forest inventory data. Application of this method for different years will provide a chance to understand the impact of forest disturbance on biomass change, and to fill in gaps in some years when the inventory data is incomplete. The mapped area covered all forested land in China, and the total forest AGB was about 12,622 Mt in 2006. The relatively fine-scaled, spatially explicit forest aboveground biomass map provides critical AGB information for forest carbon cycle studies and forest resource management.