Discovering forest height changes based on spaceborne lidar data of ICESat-1 in 2005 and ICESat-2 in 2019: a case study in the Beijing-Tianjin-Hebei region of China

Abstract

Background

The assessment of change in forest ecosystems, especially the change of canopy heights, is essential for improving global carbon estimates and understanding effects of climate change. Spaceborne lidar systems provide a unique opportunity to monitor changes in the vertical structure of forests. NASA’s Ice, Cloud and Land Elevation Satellites, ICESat-1 for the period 2003 to 2009, and ICESat-2 (available since 2018), have collected elevation data over the Earth’s surface with a time interval of 10 years. In this study, we tried to discover forest canopy changes by utilizing the global forest canopy height map of 2005 (complete global coverage with 1 km resolution) derived from ICESat-1 data and the ATL08 land and vegetation products of 2019 (sampling footprints with 17 m diameter) from ICESat-2.

Results

Our study revealed a significant increase in forest canopy heights of China’s Beijing-Tianjin-Hebei region. Evaluations of unchanging areas for data consistency of two products show that the bias values decreased significantly from line-transect-level (− 8.0 to 6.2 m) to site-level (− 1.5 to 1.1 m), while RMSE values are still relatively high (6.1 to 15.2 m, 10.2 to 12.0 m). Additionally, 58% of ATL08 data are located in ‘0 m’ pixels with an average height of 7.9 m, which are likely to reflect the ambitious tree planting programs in China.

Conclusions

Our study shows that it is possible, with proper calibrations, to use ICESat-1 and -2 products to detect forest canopy height changes in a regional context. We expect that the approach presented in this study is potentially suitable to derive a fine-scale map of global forest change.

Background

Forest vertical structures are constantly affected by natural disturbances and human activities (Vogelmann et al. 2009), causing canopy height changes. Quantitative information on the distribution of forest canopy height is essential for estimates of forest biomass and terrestrial carbon flux (Wang et al. 2016). As traditional forest inventory is very time-consuming and often limited to local scales (Bergen et al. 2009; Wulder et al. 2012), remote sensing plays an increasingly important role due to the large-scale coverage. However, passive optical remote sensing has shown limited sensitivity to forest structure above a low threshold of biomass (Blackard et al. 2008), and these sensors are primarily responsive to canopy cover, rather than vertical structure (Chopping et al. 2008). Although current endeavors to map landscape dynamics are still dominated by the use of multi-source satellite imagery data (DeVries et al. 2015), lidar provides an alternative technique to measure vertical structure changes of vegetation.

Over 50 years of research has demonstrated the utility of airborne lidar for natural resource assessment (Nelson 2013). A number of studies related to forest height changes utilized repeated airborne lidar data in local area, achieving bi-temporal or multitemporal observations (Dubayah et al. 2010; Hudak et al. 2012; Zhao et al. 2018). While airborne lidar data provides detailed vegetation structures with high accuracy, its application over large areas is constrained by the high expense associated with data acquisition (Liu et al. 2019). Spaceborne lidar possesses the capability of capturing vertical information globally and efficiently. Before the launch of Ice, Cloud and land Elevation Satellite (ICESat) -2, due to limited availability of data, it was difficult for spaceborne lidar to monitor forest change for a long period of time. Currently, combining the NASA’s ICESAT − 1 and − 2 may alleviate this problem.

The Geoscience Laser Altimeter System (GLAS) on board the ICESat-1 collected waveform lidar data from 2003 to 2009. It was the first spaceborne altimetry mission focusing on the cryosphere and also contributed to measuring global surfaces including oceans, land, and vegetation. GLAS offers an unprecedented opportunity for estimating forest canopy height and biomass from regional to global scale (Rosette et al. 2008). However, all past and planned spaceborne lidar systems collect data in the form of transects of individual footprint observations spaced at some distance from each other (Lefsky 2010). The data was not spatially continuous. To address this limitation, some researchers (Lefsky 2010; Simard et al. 2011; Wang et al. 2016) published their global forest canopy height maps with different height metrics and resolution. These maps combined spatially continuous ancillary variables to obtain complete coverage in the horizontal plane, such as Moderate Resolution Imaging Spectroradiometer (MODIS) products, climate data or terrain data. Due to map accuracy and resolution, global canopy height maps based on GLAS data present a fine-scale description of forest vertical structures at regional and global scales (Simard et al. 2011), but they are inappropriate for local applications.

With the launch of ICESat-2 in September 2018, photon counting technology has been utilized by the Advanced Topographic Laser Altimeter System (ATLAS) instrument (Neuenschwander and Pitts 2019). The advantage of the photon counting lidar is a reduced laser power requirement that allows for a smaller payload. This makes it possible to achieve a higher repetition that facilitates improved along-track spatial resolution (Popescu et al. 2018). To further improve the longitudinal spatial sampling in the mid-latitudes (approximately 60°S to 60°N) to be better than 2 km between equatorial ground tracks, ICESat-2 off-nadir points a maximum of 1.8° from the reference ground track rather than maintaining the repeat ground track coverage (Neuenschwander and Pitts 2019). Benefiting from the increased spatial sampling, ICESat-2 mission also produces geophysical products over various surface types such as land ice height (ATL06), sea ice elevation (ATL07), land/vegetation elevation (ATL08), atmosphere characteristics (ATL09), sea ice freeboard (ATL10), ocean elevation (ATL12), and inland water elevation (ATL13) (Markus et al. 2017). Among them, the ATL08 products facilitate forest assessment at global scale and promote carbon monitoring (Liu et al. 2019). Now, the second release (version 2) of ICESat-2 products is available, which has fixed a list of known issues relevant to the first release (Neuenschwander and Magruder 2019). It is anticipated that the second release will have better performance and provide more reliable data. However, few studies are concentrating on the practical applications of ALT08 products for assessing forest vertical structure.

Given the decadal gap between the two generations of spaceborne lidar, it is also an opportunity to investigate the Earth surface’s changes in a relatively long period of time. In this work, we investigated the canopy height changes in China’s capital Beijing and its surrounding areas, which is a first attempt to combine two different types and generations of spaceborne lidar data in a practical application. Unlike optical remote sensing images that mainly reflect forest cover change, or airborne lidar that mostly applied in local areas, this study utilizes spaceborne lidar to detect forest height changes within a large region. Our specific objectives are as follows: (1) to develop a universal methodology for comparing a global forest canopy height map and specific ATL08 products; (2) to evaluate the consistency of both the two products in unchanging areas; (3) to present recommendations for height change monitoring of plantation forests.

Materials and methods

Study area

The Beijing-Tianjin-Hebei Urban Agglomeration, also known as the Capital Economic Circle, is an important core area of northern China both in economic aspect and ecological environment. With approximately 218,000 km², this region consists of a northwestern ecological conservation zone, a central core functional zone, a southern functional expansion zone, and an eastern coastal development zone. Among these, the ecological conservation zone accounts for nearly half of the total area (Fig. 1). The dominant species of the mountainous part in our study area are Pinus tabuliformis, Larix gmelinii, Platycladus orientalis and Betula platyphylla. In the plain area, the dominant species are Sophora japonica, Ulmus pumila and Salix babylonica.

Fig. 1

Location of the study area (left), and the regional division (right)

For decades, tremendous efforts have been invested to improve forest coverage, forest quality, and the air quality in this important region. Numerous policies and projects were implemented, and most of them are still relevant. The “Three-North Shelter Forest Program” is the most famous one (Duan et al. 2011). The “Beijing-Tianjin Sandstorm Source Control Program”, the “Grain for Green Program”, and the “Natural Forest Protection Program” have also made great contributions in improving the quality of environmental services. For example, forest cover in Beijing increased from 35.5% in 2005 to 43.5% in 2018.

Selection of forest canopy height maps

Two aspects were mainly considered when selecting the maps: resolution and height metric. The resolution should be as high as possible, and the height metric should be consistent with ATL08 products. There are three candidate maps, published in 2010, 2011, and 2016. Lefsky (2010) produced the first map at an average resolution of 25 km² patches by combining GLAS and MODIS data. This map chose Lorey’s height as height metric, a basal area weighted height of all trees. The 2011 product of Simard et al. (2011) represented a considerable improvement by increasing the resolution to 1 km. They used climate and terrain data as ancillary variables, which are closely related to forest growth. The height metric is the waveform metric RH100 defined as the distance between the beginning and the lidar ground peak of the recorded waveform. The last product of Wang et al. (2016) went a step further by providing a spatial resolution of 500 m. They extracted the distance between the first and last peaks by using the Gaussian and wavelet methods, representing the mean height of the dominant and codominant trees as the height metric. In addition, these maps utilized different land cover products to define forest borders. All these maps have similar nominal accuracy (Table 1), but with a large disparity among them (Wang et al. 2016).

Table 1 Comparison of the three candidate maps

ICESat-2 ATL08 products contain canopy height metrics of different percentiles, and define the relative 98% height of classified canopy photon heights above the estimated terrain surface as the top of canopy. The relative 98% height rather than 100% is preferred because solar background noise may not have been fully rejected (Neuenschwander et al. 2019). Compared with the three height metrics adopted by these maps, RH100 gave the best matches. The 2011 and 2016 products have a high agreement, but the difference of mean values of the two height metrics in different forest types ranges from 5 to 7 m (Wang et al. 2016). Therefore, we selected the 2011 product of Simard et al. (2011), (henceforth, the Simard-map) as our final global map (Fig. 2).

Fig. 2

Simard et al. (2011) map of the study area

Preprocessing of ATL08 products

We selected all of the ATL08 data (National Snow and Ice Data Center, USA 2018) that covered our study area between 2019 and 05-01 and 2019-09-06 (Release 002), in which period trees have leaves for both evergreen and non-evergreen forests (total: 57 granules). Simard et al. (2011) selected GLAS data acquired from May and June of 2005. We used a longer period to obtain a higher sampling rate.

ICESat-2 have 3 pairs of laser beams, each consisting of a strong and a weak energy beam. Strong beams are recommended for lower reflectance objects (e.g., vegetation), but weak beams still provide valuable data especially when the background noise rate is small in night acquisitions (Neuenschwander and Magruder 2019). Figure 3 shows the relative frequency distributions of height values of strong and week beams. Since strong beams are more reliable, we deleted all the weak beams (Arendt et al. 2019).

Fig. 3

Relative frequency distributions of height values of strong beams (left) and weak beams (right)

By labeling ground photons and canopy top photons sequentially, the ATL08 algorithm extracts terrain and canopy height from ATL03 (geolocated photons product) photon clouds. If there are fewer than 50 labeled signal photons in the 100 m step, no calculations will be made and invalid values will be reported in the dataset (Popescu et al. 2018). Sometimes, mislabeling occurred causing outliers (Neuenschwander and Magruder 2016). Forest heights of the Simard-map range from 5 to 39 m in our study area and pixels between 36 and 39 m accounting for only 0.14% (119/84042). Extreme values were discarded because they are less reliable. Considering the low vegetation in barren areas, we decided that 2 m (instead of 5 m) to 35 m to be reasonable for forest canopy height assessment in our study area. The invalid values and outliers were filtered out.

Data comparison method

It is important to note that the canopy height map and ATL08 products represent different types of data with different scales. In fact, the map has a constant resolution of 30 arc sec (approximately 1 km), whereas ATL08 products provide canopy information at 100 m, in fixed-length steps along the ground track. Thus, each pixel of the map represents the maximum canopy height within a 1 km × 1 km area, while each point of ATL08 product represents a 100 m × 17 m area (the diameter of ICESat-2 footprints is 17 m). ATL08 products should be transformed into a 1 km × 1 km scale. Geographical matching of ATL08 points and the map pixels was thus completed first.

Regarding ICESat-2 track lines which allow systematic sampling on the earth’s surface, each pixel intersected by the track lines has about 10 sampling points ideally. We only preserved the sample point with a maximum value in one pixel, which represents the maximum canopy height in 2019. For a small number of pixels, few points remain after filtering. These pixels and points were abandoned, because the maxima are unrepresentative with too few sample points. We defined the threshold as 3 points, because it is a proper balance between data volume and representativeness (Fig. 4). Thus, we obtained the maximums in 1 km × 1 km scale. Finally, the difference between the pixel values and the maxima indicates the forest canopy height changes at the 1 km × 1 km scale between 2005 and 2019 (Fig. 5). Unchanging areas (like protected natural forests) and changing areas with three major patterns (tree growth, forest disturbance and new plantation) are classified for further evaluations.

Fig. 4

Percentage of preserved maxima under different thresholds of point quantity

Fig. 5

Diagram of the comparison methodology (the distance between the track lines of ATL08 products was shortened for improved interpretation). Strong beams were selected and filtered from the original ATL08 products, and geographically matched with the map. The change of canopy height was derived from the pixel value and the maximum of ATL08 points in one pixel. Protected natural forests are regarded as unchanging areas. Areas where change can be expected are divided into three patterns

Evaluation method

Our evaluation focused on whether this processing method can correctly reflect the forest canopy height changes, instead of the accuracy of each product. For new plantations and disturbed forests, we ought to detect the difference; for protected natural forests, no difference (or only a minor change) is expected. We assume that the canopy heights in protected natural forests are constant within a local area of 1 km × 1 km resolution. Thus, we regard protected natural forests as unchanging, and all other areas as “changing” where change can be expected.

Evaluating changing areas at the 1 km × 1 km scale is a challenging task. We do not know the exact situation in 2005. Thus, we decided to evaluate the consistency of the unchanging areas instead of the difference in the changing areas. Unchanging areas are usually associated with high terrain slope, where the combination of the vegetation signal and the terrain signal results in a high uncertainty for GLAS data of ICESat (Lefsky et al. 2007; Hilbert and Schmullius 2012). Benefitting from the small footprint and high laser repetition rate, ICESat-2 laser technology overcomes this problem. If our processing method could get acceptable results in topographic relief areas, we can expect even better performance in relatively flat areas.

We selected three unchanging sites, as shown in Table 2 and Fig. 6. They are all protected natural forests with continuous, dense, and tall trees. At the ICESat-2 line level, the average bias and root mean square error (RMSE) between the pixel values and the maximums of each intersected line were calculated, noted as line-transect-level. Further, the bias and RMSE were averaged for each site, noted as site-level. By changing the sampling point threshold from 1 to 10, we also tested its impact on bias and RMSE at site-level.

Table 2 Evaluation sites information

Fig. 6

Maps of the three evaluation sites (Intersected transects are shown in yellow lines)

Results

Comparison results

Before comparing the data sets, forest canopy heights derived from every filtered ICESat-2 step (total: 170,538 points) are shown in Fig. 7a (height: 2–35 m). The northern region preserves more forests than the southern region, while the central urban area is surrounded by plenty of tall forests. Figure 7b presents the changes of forest canopy height between 2005 and 2019 at a 1 km × 1 km scale. Since Simard-map (2011) has the limitation that areas with tree heights less than 5 m are regarded as 0 m (noted as ‘0 m’ pixels), we ignored the maxima of ATL08 products less than 5 m resulting in a total of 20,045 maxima preserved. In general, forest height increases exceed forest height decreases (74% vs. 26%; 14,900 pixels vs. 5,145 pixels). The increasing areas are mostly distributed in the northern and middle region, and the southern region changed slightly except for the western edge.

Fig. 7

Forest canopy heights derived from filtered ICESat-2 points (a), changes of forest canopy height between 2005 and 2019 (b), and height changes of ‘0 m’ pixels (c)

An interesting discovery is the fact that more than half (58%; 98,076/170,538) of ATL08 points lie in ‘0 m’ pixels with an average height change of 7.9 m. If considering maxima at the 1 km × 1 km scale, the average change is 13.7 m. Figure 7c demonstrates the changes in ‘0 m’ pixels, which are more likely to reflect the tree plantation distribution. Height changes mainly occur in the ecological conservation zone, where various tree planting programs have been carried out.

Evaluation results

In Fig. 7, the three protected forests, defined as unchanging area, are also marked as three red boxes. Within the boxes, the color is generally dark green in Fig. 7a and light yellow in Fig. 7b. Table 3 lists the bias and RMSE for the canopy height map compared with the ATL08 products for both line-transect-level and site-level. For all the lines, the bias values range from − 8.0 to 6.2 m, and the RMSE values range from 6.4 to 15.1 m. It should be mentioned that the accuracy evaluation of forest height products in mountain areas is generally missing (Simard et al. 2011; Wang et al. 2016). The major reasons include the limited ground truth data and potential high errors. Thus, the high bias and RMSE errors at line-transect-level are expected and understandable in mountain areas. At the site-level, the bias values range from − 1.5 to 1.1 m, and the RMSE values range from 10.2 to 12.0 m. The RMSE errors are still relatively high, but the bias decreased significantly. Figure 8 shows all the 769 pixels values and maximums values.

Table 3 Evaluation results of the three unchanging sites

Fig. 8

Comparison of canopy heights between 2019 and 2005

Table 4 shows the bias and RMSE varying with the sampling point threshold. Generally, the RMSEs almost remain stable while the biases fluctuate slightly with increasing threshold value (Table 4). The reason why the threshold setting is not reflected in the results is that we chose continuous and dense protected forests as evaluation sites. These areas are homogenous such that the difference in the representativeness of one point as against ten points is negligible. However, we still recommend setting a reasonable threshold value to get a proper balance of data volume and representativeness. For most regions, it is expected that the increase of points in one pixel can reduce the error caused by ‘sampling’.

Table 4 Evaluation results under different threshold values

Greening examples

Figure 9a shows the first comparison example (115.88°E, 40.48°N) of a line transect of approximately 3 km, which represents a typical site of China’s Grain for Green Program. The left two panels depict the location of the ICESat-2 track in Google Earth. This site used to be farmland in 2007, and now is green forest. Consistent with the images, intersected pixels of the canopy height map display ‘0 m’ here in 2005, while ICESat-2 ATL08 products provide more detailed information about the forest canopy height in 100 m segments. In the upper part of the track, a small grass area (light green) in this line has no height value. By comparing the two-generation data, we can discover the new plantation areas, and then obtain details from ATL08 products.

Fig. 9

Two quantitative greening examples. a The left two panels depict the location of the ICESat-2 transect in Google Earth for context (acquired on 2007-09-21 and 2019-10-02), the right panel shows the comparison of the canopy height map data and ATL08 product (ATL08_20190815051930_07350406_002_01). b The Google Earth context data were acquired on 2008-02-20 and 2019-07-30, and the product is ATL08_20190508095629_06130306_002_01

The second line example (117.61°E, 42.17°N) is also an approximately 3 km transect, but in a mountain region (shown in Fig. 9b). The center of the left panel image is covered by several patches of forests, while other areas are bare earth. We acknowledge that this picture taken in winter is not really suitable for some non-evergreen trees which may be masked in the brown soil. We use it nevertheless because no other images in Google Earth could meet both clarity and time requirement, and the canopy height map also displays the same forest feature. Over the past decade, the central forest area increased and is now covering almost the whole region. Referring to ALT08 products, the canopy height map obviously overestimates the height of the central forest, which is an inevitable error caused by the global map production.

Discussion

The goal of this paper is to explore the possibility of comparing the two generations of spaceborne lidar data and to attempt to utilize the ICESat-2 data in a practical application. In this study, we used a 1 km × 1 km scale to avoid two difficulties: small overlap probability (footprint locations of ICESat-1 and -2) and regional inconsistency (70 m diameter for ICESat-1 and 100 m × 17 m for ICESat-2). We obtained the canopy height distribution in 2019, the canopy height changes between 2005 and 2019, and the ‘0 m’ pixel changes of the Beijing-Tianjin-Hebei Urban Agglomeration, China. A recent study shows that China and India are leading in the “greening of the world” through land-use management (Chen et al. 2019). Our study also confirmed this finding from the aspect of forest canopy height change and distribution. Limited by the quantity of valid transect data (meet season, cloud, beam power and other conditions), our results now seem a little sparse that only reflect overall status. But we believe this problem will be solved by the increase of satellite on-orbit time and eventually depict fine maps.

Although we took measures to reduce errors, there are several factors that contributed to the uncertainty of the comparison results. First and foremost, a proper comparing method may be the most challenging procedure, as this part induces extra errors and could even cause unrealistic results. We used one of these global canopy height maps, which has continuous coverage, as our base map of 2005. Comparing with ICESat-2 data confronted two choices — direct comparison or map to map comparison. The first method is easier to operate, but has the discrepancy problem of data type and also results in discontinuous transect lines. The second method is relatively complex that needs extrapolating to continuous map like previous canopy height maps before comparing. In this way, weaknesses of the first method can be avoided. But this process induces not only the comparison error, but also the map production error. Furthermore, producing a continuous canopy height map based ICESat-2 data is not the major purpose of this paper. Therefore, regarding ICESat-2 track lines as systematic sampling, we preserved the maximum point in one pixel which has 3 or more points and compared the maximums with pixel values. It is acknowledged that several sampling points were not completely representative, but we consider it was acceptable. Benefiting from the off-pointing feature, discontinuous track lines will become more and more dense that could also reflect relatively detailed local canopy height information.

Another factor that impacts the accuracy is the difference of land cover products employed by the canopy height map and the ICESat-2 ATL08 product. The canopy height map used GlobCover map (Hagolle et al. 2005) to define forest area, while ICESat-2 utilized Landsat Vegetation Continuous Fields (VCF) product (Sexton et al. 2013). It should be noted that the canopy existence is determined in a 10-km segment in the ATL08 algorithm. If the average VCF canopy cover along each 10 km segment is > 5%, canopy is presumed along the segment and processed accordingly (Neuenschwander and Pitts 2019). It is hard to determine the exact impact of this difference, since land cover types also changed in the past decade years. But we suspect different products and processing methods inevitably caused boundary discrepancy in edge and local areas of forest, resulting in abnormal height changes. In cases where the GlobCover map defined non-forest, ATL08 algorithm may detect more delicate forest areas. Thus, ‘0 m’ pixel changes may not fully represent tree planting program efforts. We can only conclude that these changes are more likely to reflect the tree plantation status.

Besides the uncertainty induced by the comparison process, possible errors should be considered. Areas with very low or high canopy heights usually involve greater uncertainty (Bolton et al. 2014). Moreover, in mountain areas, high topographic relief caused greater uncertainties in the global map. Due to the nature of photon counting lidar and the ATL08 algorithm, there could be instances where the noise photons are misidentified as either canopy or terrain photons, which subsequently leads to incorrect estimates of canopy heights, especially for extremely dense or sparse vegetation areas (Neuenschwander and Pitts 2019). In this study, we selected three protected forests as unchanging sites to evaluate the agreement among both. The results show that bias values decreased significantly from line-transect-level to site-level, while RMSE values are still relatively high. It is expected that our comparison results perform better in low topographic relief areas with median canopy cover.

Conclusions

This study presents an investigation of the feasibility of comparing two generations of spaceborne lidar data to discover forest canopy height changes. A global forest canopy height map and available ATL08 products were utilized, representing the forest status in 2005 and 2019. Our results indicate a remarkable increase in forest canopy heights which is mostly caused by China’s ambitious planting programs. Overall, our comparison methodology provides a new opportunity to monitor the forest change on a global scale, especially regarding new plantations. We are convinced that by combining ICESat-1 and -2 data in vegetation change studies, more accurate results can be obtained, especially when more sophisticated procedures are developed. Further research is required to improve sampling coverage of ICESat-2 data with weak power beams. With the increasing availability of ICESat-2 data and field data, we expect that our method will improve global forest change mapping.