Retrieval of forest canopy height in a mountainous region with ICESat-2 ATLAS

Shiyun Pang, Guiying Li,＊, Xiandie Jiang, Yaoliang Chen, Yagang Lu,Dengsheng Lu

a State Key Laboratory for Subtropical Mountain Ecology of the Ministry of Science and Technology and Fujian Province,Fujian Normal University,Fuzhou,350007,China

b Institute of Geography, Fujian Normal University, Fuzhou, 350007, China

c Institute of East China Inventory and Planning, National Forestry and Grassland Administration, Hangzhou, 310019, China

Keywords:Airborne LiDAR Canopy height ICESat–2 ATLAS Mountainous region Segment filtering

ABSTRACT Forest canopy height is one of the important forest parameters for accurately assessing forest biomass or carbon sequestration. ICESat-2 ATLAS provides the potential for retrieval of forest canopy height at global or regional scale, but the current canopy height product (ATL08) has coarse resolution and high uncertainty compared to airborne LiDAR-derived canopy height (hereafter ALCH) in mountainous regions, and is not ready for such applications as biomass modeling at finer scale. The objective of this research was to explore the approach to accurately retrieve canopy height from ATLAS data by incorporating an airborne-derived digital terrain model(DTM) and a data-filtering strategy. By linking ATLAS ATL03 with ATL08 products, the geospatial locations,types,and(absolute)heights of photons were obtained,and canopy heights at different lengths(from 20 to 200 m at 20-m intervals) of segments along a track were computed with the aid of airborne LiDAR DTM. Based on the relationship between the numbers of canopy photons within the segments and accuracy of ATLAS mean canopy height compared to ALCH,a filtering method for excluding a certain portion of unreliable segments was proposed.This method was further applied to different ATLAS ground tracks for retrieval of canopy heights and the results were evaluated using corresponding ALCH. The results show that the incorporation of high-precision DTM and ATLAS products can considerably improve the retrieval accuracy of forest canopy height in mountainous regions.Using the proposed filtering approach, the correlation coefficients (r) between ATLAS canopy height and corresponding ALCH were 0.61–0.91, 0.65–0.92, 0.68–0.94 for segment lengths of 20, 60, and 100 m, respectively;RMSE were 1.90–4.35, 1.55–3.63, and 1.34–3.23 m for the same segment lengths. The results indicate the necessity of using high-precision DTM and using the proposed filtering method to retrieve accurate canopy height from ICESat-2 ATLAS in mountainous regions with dense forest cover and complex terrain conditions.

1. Introduction

Forests, as the largest carbon pool in the terrestrial ecosystems, play important roles in protection of ecological environments, global carbon cycling,and reduction of climate change(Pourshamsi et al.,2021;Wang et al., 2021). Canopy height is one of the most important forest parameters for modeling forest biomass or carbon stock (Tang et al., 2019;Narine et al., 2020). Since the ICESat (Ice, Cloud, and Land Elevation Satellite) GLAS (Geoscience Laser Altimeter System) data in 2003 and recent ICESat-2 ATLAS (The Advanced Topographic Laser Altimeter System)and GEDI(Global Ecosystem Dynamics Investigation)LiDAR are available,many studies in the last two decades have explored methods to retrieve canopy height at national or global scale (Potapov et al., 2020;Liu et al.,2022).However,how to quickly and accurately extract canopy parameters from satellite LiDAR data is still challenging, especially in mountainous regions.

Traditionally,forest canopy height is obtained from field surveys,but the high cost for fieldwork, uncertainty in tree height measurement in dense forests (especially for tall trees in mountainous areas), and the difficulty in updating data constrain its extensive application at large scales. Airborne or UAV (Unmanned Aerial Vehicle) LiDAR can provide 3-dimensional features, and have advantages over field measurements,thus have gained attention in recent years for retrieval of tree attributes or forest structure parameters (Lu et al., 2020; Wang et al., 2020).However,its large data volume and high cost for data acquisition limit its applications in small areas (Hayashi et al., 2013; Liu et al., 2019).Compared to airborne LiDAR, spaceborne LiDAR provides global coverage at relatively lower cost, has potential for retrieving forest height,and can be used for biomass or carbon stock modeling at regional or global scale(Enβle et al.,2014;Pang et al.,2019;Sun et al.,2020).The first LiDAR satellite,ICESat GLAS,successfully accomplished its mission in 2010, providing global multi-year topography and vegetation data during 2003–2009, from which various vegetation products such as canopy heights were derived and used for global or regional forest biomass modeling (Lefsky, 2010; Healey et al., 2015; Santoro et al.,2021). Compared to ICESat GLAS, the improved features of ICESat-2 ATLAS have many advantages in canopy height retrieval, especially in complex terrain conditions (Narine et al., 2019a, 2019b; Neuenschwander et al., 2019a; Duncanson et al., 2020), and may provide new opportunities for extracting accurate canopy parameters.

Previous studies using ICESat-2 ATLAS data for retrieval of forest canopy height were mainly for plantations or forests in flat areas(Li et al.,2020; Lin et al., 2020; Zhu et al., 2020; Jiang et al., 2021). Undulating terrain conditions with dense forest cover can seriously affect the number of photons reaching the ground and the precision of photon classification,thus further influence the extraction of digital terrain model(DTM)and canopy height model (CHM) data (Neuenschwander et al., 2019b;Dong et al., 2021). Some previous studies indicated that the root mean square error (RMSE) between ICESat-2 DTM and airborne LiDAR DTM can be less than 1 m in flat regions(Xing et al.,2020),while the RMSE of forest canopy height can be as high as 4.55–5.56 m in hilly areas and 7.02–8.92 m in mountainous areas,and the major errors are from the low accuracy of DTM (Zhu et al., 2018). As terrain slope increases, the accuracy of ICESat-2 DTM becomes worse (Wang et al., 2019). Similar research in Canada also indicated that correlation coefficient (r) and RMSE between the canopy heights from ATLAS and airborne LiDAR were 0.84 and 2.9 m, respectively, in flat areas, but r decreased to 0.64 and RMSE increased to 3.5 m in complex forest stand structures (Queinnec et al., 2021). Forest cover also greatly impacts the retrieval accuracy of DTM from ICESat-2.As forest cover increases,the accuracy of DTM and canopy height derived from ATLAS decreases. The RMSE between ICESat-2 DTM and airborne LiDAR DTM can be 1.98 m in areas with low canopy cover (25%–40%) and 5.54 m in areas with high canopy cover(83%–84%) (Huang, 2021). Dong et al. (2021) found that the RMSE between the average canopy heights from ICESat-2 and airborne LiDAR in temperate forests and tropical rain forests were 2.55 and 4.26 m,respectively. These studies indicated the difficulty in retrieval of forest canopy height from ATLAS in complex mountainous regions or dense forest cover areas. However, many forests are located in mountainous areas with high canopy cover. To take full advantage of ICESat-2 dense LiDAR footprints,our task was to develop a suitable approach to extract canopy height from ATLAS in mountainous regions.

The accuracy of canopy height retrieved from ATLAS is also affected by the size of the statistical unit, i.e., the length of segment along the track direction(Silva et al.,2021).The current ATLAS vegetation height product ATL08 contains statistical parameters of canopy and terrain based on a fixed 100-m length.The rigid 100-m distance is too coarse and not practical for many applications, especially when finer-resolution optical or microwave data are combined for wall-to-wall canopy height interpolation or biomass modeling. Evaluation of canopy height inversion from ATLAS data in terms of segment length showed that the best lengths vary, depending on the study sites, especially the topographical conditions. Chen et al. (2019) found the best length to be 50 m in pit-and-mound topography,while Huang(2021)found 20 m was suitable for most of nine study areas with various topographic conditions (from flatland to low hill to hill), forest types (boreal, temperate, subtropical,and tropical),and vegetation covers(ranging between 1%and 87%).At a coarse level,Li et al.(2020)identified 250 m as the optimal length for the Inner Mongolia Autonomous Region of China. In recent years, most studies for retrieval of canopy height from ATLAS data used 30 m for linking it with Landsat or Sentinel-2 data for wall-to-wall mapping of forest canopy attributes (Liu et al., 2019; Li et al., 2020; Jiang et al.,2021).To date,no general conclusion about the best segment length has been obtained because of compound impacts of forest type, vegetation cover,and terrain condition.

Retrieval of canopy height from ATLAS data involves photon or segment screening, removal of unreliable photons and segments to ensure high confidence in the derived canopy heights. In general, six approaches have been developed,including(1)set slope thresholds and keep only segments with slopes falling in the defined thresholds, for example,less than 25°(Zhu et al.,2020);(2)use descriptive information in the metadata, for example, exclude the photons flagged with clouds(Neuenschwander et al., 2020) or segments with signal-to-noise ratio greater than 6(Zhu et al.,2020);(3)set a canopy height range based on prior knowledge about the forest height in a study area, keeping all segments within the defined range, for instance, 2–35 m based on field measurements(Sun et al.,2020)and 2–60 m according to reference tree height data(Li et al.,2020);(4)use only ATLAS's strong beams to retrieve forest canopy height, considering weak beams are prone to impact by forest canopy cover (Sun et al., 2020; Xing et al., 2020); (5) specify the maximum difference between a photon's absolute height and reference height,and exclude the photons when the difference is greater than the specified value,for example,20 m(Huang et al.,2020);(6)based on the number of photons/segments falling within a defined pixel size,exclude the pixels containing, for example, less than three ATL08 segments within 1000 m×1000 m(Sun et al.,2020),or keep the pixels that have 40–170 photons within 30 m×30 m(Zhu et al.,2020),or remove pixels containing less than four ground photons within 30 m×30 m(Lin et al.,2020).These photon or segment-filtering approaches are only effective in a specific area at a specific spatial scale where the approach is developed,but not suitable for regions having complex topography and dense forest cover (Dong et al., 2021). Therefore, a filtering method to ensure both quantity and quality of the photons/segments for retrieval of canopy height in mountainous regions is needed.

To evaluate the accuracy of ICESat-2 ATLAS canopy height products in a mountainous area,we created a graph profiling the photons from the ATL03 product and airborne LiDAR as well as canopy height from ATL08 and airborne LiDAR of a 1-km-long segment along the ICESat-2 track(Fig. 1). It clearly shows the discrepancy between ATLAS products and airborne LiDAR data,especially for attributes of the canopy elevation and ground elevation.The differences of absolute mean canopy heights from ATL08 and airborne LiDAR are closely related to topographic conditions and forest covers, which greatly impact the quantity and quality of photons hitting the ground and canopy. As phase I in Fig. 1 shows, the dense forest cover constrains the photons going through the canopy,resulting in few photons reaching the ground,while most photons remain on top of the forest canopy;in relatively sparse forest canopy(phases III and IV in Fig. 1), many ground photons are obtained, providing the capability of accurately extracting topographic features.In complex and undulating topography(phase II in Fig.1),some photons appeared over the canopy,resulting in spurious canopy results.This example indicates the difficulty in retrieving accurate canopy height from ATLAS data alone in mountainous regions with dense forest cover,requiring incorporation of an extra data source to improve the retrieval accuracy.

The overall goal of this research was to explore an approach to retrieve accurate canopy height in a complex topographic region with dense forest cover through a comparative analysis of using ICESat-2 ATLAS at different segment lengths, including development of an optimal filtering method to exclude low-quality photons/segments.Specifically, the objectives are to (1) evaluate ATL08 products using airborne LiDAR-derived DTM and digital surface model(DSM)data;(2)extract canopy heights by incorporating airborne LiDAR DTM and ATLAS data at different segment lengths to improve accuracy of retrieved canopy height; (3) determine the criteria of segment filtering at different segment lengths,ensuring that as many high-quality segments as possible are retained for future use; (4) validate the proposed segment filtering method by applying it to different ground tracks.This research provides new insights for retrieval of accurate canopy height at flexible segment lengths in mountainous regions when precision DTM is available. Actually, airborne LiDAR data are available in many areas, but due to the expensive data acquisition,it is unfeasible to obtain it multiple times;in turn,airborne LiDAR is rarely used for forest monitoring.However,DTM does not change over time, in general. Therefore, incorporation of satellite-derived DSM and LiDAR-derived DTM becomes possible to timely update CHM.This kind of canopy height data is especially critical for accurate forest biomass or carbon stock modeling in a large area.

Fig. 1. The profile of ICESat-2 ATLAS and corresponding airborne LiDAR data along the track. Gray and black dots represent the airborne LiDAR points; green, red,and yellow dots represent the ALT03 photons obtained by linking them to the ATL08 product; blue and purple dots represent the absolute mean canopy height from ALT08 and airborne LiDAR with segment length of 100 m.Phases I,II,III,and Ⅳrepresent dense forest cover,undulating terrain,relatively short canopy height,and relatively sparse canopy. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Fig.2. The locations of study sites in China and ICESat-2 ATLAS ground tracks.(a)Two sites in Anhui Province;(b)study site in Jinzhai County with overlaid satellite ground tracks; (c) study site in Huangshan District with overlaid satellite ground tracks.

2. Materials and methods

2.1. Study area

Two sites in the northern subtropical mountainous region of China—Jinzhai County and Huangshan District in Anhui Province — were selected as the study areas (Fig. 2). Jinzhai County is located at the western border of Anhui,the hinterland of Dabie Mountains,adjacent to Henan and Hubei provinces to the west.With a total area of 3814 km2,Jinzhai is the largest county in Anhui Province. The topography is characterized by rough terrain with an average elevation of 500 m and an average gradient of 21%. The climate is typical northern subtropical humid monsoon,with annual average temperature of 12.3°C and annual precipitation of 1600 mm.The forestlands cover 2948 km2of the county,and typical forest types are mixed deciduous and evergreen broadleaf(Gu et al., 2014). The high elevation variation (from 30 to 1663 m) in this county results in obviously vertical zonal distribution of vegetation types,comprised of evergreen broadleaf forest, mixed forest of evergreen and deciduous broadleaf, and deciduous broadleaf forest from low to high elevation.The main tree species include Pinus massoniana,Cunninghamia lanceolate, Cyclobalanopsis glauca, Quercus variabilis, Quercus acutissima Carruth., Castanea mollissima, Liquidambar formosana Hance, Carya cathayensis, Phyllostachys heterocycla (Carr.) Mitford cv. Pubescens (Fu,2018). Based on 107 forest plots surveyed in 2019, the average canopy height of the upper story was about 10 m (2.9–18.7 m), the average canopy cover was 0.68(0.25–0.90),the average understory shrub cover was 0.3 (0.03–0.80), depending on the forest types and ages and elevation.

Huangshan District is located in the south of Anhui Province,including the famous Huangshan scenic area.The total area is 1775 km2with forestland of approximately 1234 km2(Cui, 2015). The region is characterized by raised peaks and plunging ravines,the elevation ranges from 64 to 1806 m with 23%of total area having slopes of over 25°.This region has a humid subtropical moist monsoon climate with annual average temperature of 15.4°C and precipitation of 1500–1600 mm,with most occurring in spring and summer. The unique topographical and climatic merits created a remarkable habitat that facilitates great biodiversity. The typical forest type is evergreen broadleaf forests.However, due to the peculiar terrain, the vegetation on Huangshan Mountain also has a distinctive vertical distribution: the plants at the high, middle, and low elevations belong to the frigid, temperate, and subtropical biomes, respectively. The main tree species include Pinus massoniana, Pinus taiwanensis, Cunninghamia lanceolata, Phoebe sheareri,Cinnamomum camphora, Lithocarpus glaber, Castanopsis sclerophylla,Daphniphyllum macropodum, Cyclobalanopsis glauca, and Liquidambar formosana Hance.

2.2. Datasets used in research

The airborne LiDAR point cloud data were acquired using the RIEGL VQ-1560i-DW LiDAR scanning system mounted on an airplane in October 2018 for Huangshan District and June 2019 for Jinzhai County.The LiDAR system uses a near-infrared spectral channel,and LiDAR point cloud density is 4 points·m-2in Huangshan and 2 points·m-2in Jinzhai.Before their delivery,the LiDAR data were de-noised(the isolated points and noise were deleted), and LiDAR points were labeled as ground and non-ground by a data provider using the Improved Progressive Triangular Irregular Network(TIN)Densification algorithm,then projecting to CGCS 2000 3-Degree Gauss Kruger zone 39 in LiDAR 360. The ground points and the first return non-ground points were respectively used to generate DTM and DSM with cell size of 1 m using the TIN interpolation method.CHM was calculated by subtracting the DTM from the DSM and refined by removing abnormal values such as negative values and very large values by median filtering. Airborne DTM, DSM, and CHM have a normal height system.

The ICESat-2 ATLAS ATL03(the global geolocated photon cloud)and ATL08 (Land and Vegetation Height) products were downloaded from National Snow & Ice Data Center (https://nsidc.org/data/ICESat-2).Specifically, we collected ICESat-2 data from ground tracks that passed over Jinzhai between June and September 2019 and Huangshan in October 2019,close to the dates of the airborne LiDAR data acquisitions.The ATL03 product contains the parameters of latitude, longitude, and emission time of individual photons received, as well as their height above the WGS84 reference ellipsoid(Neumann et al.,2021).It serves as the input data for each of the higher-level data products. The ATL08 product is developed from ATL03 using DRAGANN (Differential,Regressive, and Gaussian Adaptive Nearest Neighbor), an algorithm specifically designed for extracting terrain and canopy heights. It provides estimates of terrain and canopy height and canopy cover of a segment at a 100-m step size in the along-track direction (Neuenschwander et al.,2021).During the processing of ATL08,the photons used for calculating the parameters of terrain and canopy heights were categorized as noise(value=0),ground(value=1),canopy(value=2),and top of canopy(value=3),and were indexed back to the ATL03 product.Based on these features, the photons in ATL08 were traced back to ATL03, and their geographic coordinate information and absolute elevation were extracted from ATL03. Previous research indicated that the signals from ATLAS's weak beams are too weak to determine terrain and canopy height, especially in densely vegetated areas (Neuenschwander et al.,2020;Xing et al.,2020);thus,only strong beams were used in this research.

2.3. Methodology

The framework of this study is illustrated in Fig. 3, including the following steps:(1)evaluation of ATL08 products using airborne LiDARderived metrics,in terms of photon types,DSM,DTM,and forest canopy heights; (2) calculation of canopy heights at different segment sizes by combining LiDAR DTM with the absolute height of canopy photons from ATL03 and ATL08 products; (3) analyses of the scale effect on canopy height retrieval by comparing it with the airborne LiDAR-derived counterpart, and determination of the optimal filtering threshold for a specified segment size;(4)validation of the proposed method by applying it to different ATLAS ground tracks.

2.3.1. Evaluation of the ATL08 product

Taking the ATLAS track GT1L on June 10,2019 in Jinzhai County as a case study,the ATL08 product quality was evaluated by two approaches:(1) We compared the profiles of ATLAS photons with airborne LiDAR clouds, especially at four locations with different slope steepness and forest covers.Because ATLAS elevation has a geodetic height system,and airborne LiDAR data has a normal height system, ATLAS elevation was first calibrated to the normal height system using the Earth Gravitational Model 2008 (EGM2008) (Berninger and Siegert, 2020). Here slope and forest cover were computed from airborne LiDAR DTM and CHM. The forest cover was defined as the ratio of number of cells with CHM higher than 5 m divided by the total number of cells within a segment (Dong et al., 2021). (2) We examined the correlation between corresponding canopy height metrics derived from the ATL08 product and airborne LiDAR. Considering that ATL08 land and vegetation height metrics are generated at a segment length of 100 m and ATLAS scanning width is 14 m,the corresponding airborne LiDAR metrics from LiDAR CHM are also calculated within 100 m × 14 m, which is spatially coincident with the ATL08 segment. These metrics include maximum, minimum, median,and mean canopy height (Hmax, Hmin, Hmed, and Hmean) and percentiles at 25, 50, 60, 70, 80, 90, and 98 (RH25, RH50, RH60, RH70,RH80, RH90, and RH98). The correlation coefficient (r) between the corresponding metrics from the two datasets were calculated and used to evaluate the quality of ATL08 products.

Fig. 3. The framework for evaluation of ICESat-2 ATL08 product and development of an approach to retrieve canopy height. AL-CHM, AL-DSM, and AL-DTM represent canopy height model, digital surface model, and digital terrain model, respectively, from airborne LiDAR data; ATL03 and ATL08 represent different levels of ATLAS products: ATL03, global geolocated photon cloud; ATL08, land and vegetation height.

2.3.2. Incorporation of airborne LiDAR DTM and ATLAS photon data for retrieval of canopy height, and determination of segment-filtering criteria at different segment sizes

A comparative analysis of ATL08 products with the counterparts of airborne LiDAR showed large differences between the two datasets,and these discrepancies are mainly due to the uncertainty of the DTM used to calculate ATL08 parameters, as shown in Fig. 1. Thus, we introduced airborne LiDAR DTM to replace the original DTM used for vegetation height calculation in the ATL08 product. In other words, we subtracted airborne LiDAR DTM from ATLAS elevation (i.e. absolute height) of canopy photons, and obtained the relative heights of individual canopy photons, from which the mean canopy height of a segment at the predefined size was calculated by averaging total relative heights over the number of canopy photons falling within the segment. To examine the effectiveness of this approach,r and RMSE between the resulting ATLAS mean canopy height and airborne LiDAR-derived mean canopy height at segment length of 100 m were calculated and compared.

The reliability of such retrieved canopy height is affected by the number of canopy photons falling within the segment.The more photons contained within a segment, the higher the accuracy of height metrics calculated from the photons; however, the fewer the qualified segment retained in a given study area.Therefore,a balance between the accuracy and data volume of high-quality segments is needed and can be obtained by identifying an optimal threshold, which is the minimum number of canopy photons within a predefined segment size required for accurately calculating height metrics,while keeping as many high-quality segments as possible.

Again taking the ATLAS track GT1L on June 10, 2019 in Jinzhai County as an example, we set the segment length at 20 m to depict the procedure of determining the threshold for selection of canopy photons:(1) compute the mean canopy height from airborne LiDAR for each corresponding segment; (2) start with a threshold value of 1, select the segments that contain canopy photons equal to or more than the threshold, count the number of such segments, and compute the mean canopy height of each ATLAS segment;(3)take the mean canopy height from airborne LiDAR as reference, calculate r and RMSE between two mean canopy heights;(4)iterate this process by increasing the threshold by 1 each time until the threshold reaches the maximum number of photons; (5) construct a graph using the threshold as the x-axis, RMSE and r as one y-axis, and the number of segments that contain canopy photons greater than the threshold (data volume) as another y-axis(Fig. 4a). Fig. 4a shows that r between two mean canopy heights of ATLAS and airborne LiDAR becomes larger as the threshold value increases,and RMSE inverses.The number of segments appears to decline with slow-fast-slow trends. Based on this characteristic, we designed a measurement called “data volume declining ratio” to evaluate the declining rate of qualified segments as the threshold value increases.Suppose the threshold values of canopy photons are X1,X2,X3,…,Xmaxin ascending order and the corresponding numbers of qualified segments are Y1,Y2,Y3,…,Ymax,then Y1–Y2,Y2–Y3,…,Ymax-1-Ymaxare reduced segment numbers when the threshold value increases from X1to X2,X2to X3,…,Xmax-1to Xmax.The data volume declining ratio is defined as the reduced number divided by the maximum reduced number(Max(Y1–Y2,Y2–Y3, Y3–Y4, …, Ymax-1- Ymax)) (Fig. 4b). In order to obtain a high r value between two mean canopy heights of ATLAS and airborne LiDAR while keeping a sufficient amount of data (high-quality segments), we can select a threshold around the peak to prevent the sharp decline of data.The detailed process to determine the threshold can be described as follows: the threshold should be between when the data volume declining ratio reaches 0.5 the first time and the second time,as shown by the red line in Fig.4b.Within this range,we selected three threshold values for comparative analysis: (1) the starting threshold (S) when the data volume declining ratio reaches 0.5 the first time; (2) the median value of the threshold (M); and (3) the retention rate (the ratio of the number of qualified segments to total number of segments) close to but not less than 20% (E). Finally, the same approach is used to determine optimal thresholds for different segment sizes,from 20 to 200 m with a step of 20 m.

Fig.4. The approach to select thresholds(example segment length is 20 m):(a)the relationships of thresholds with correlation coefficient (r), root mean squared error(RMSE)based on the mean canopy height of ATLAS and airborne LiDAR, and number of segments containing photons equal to or more than the threshold; (b)the declining ratio based on number of segments when threshold value incrementally increases by 1.The red line in(b)indicates the range of the potential thresholds;S,M,and E represent three methods for threshold selection based on the declining ratio reaching 0.5 the first time, the median value, and the retention rate of higher than 20%. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

2.3.3. Validation of proposed ATLAS canopy height retrieval approach by applying it to other ground tracks

In order to examine the robustness of the proposed approach for retrieval of canopy height through integration of airborne LiDAR DTM and ATLAS products and segment filtering, the proposed approach was applied to different ATLAS tracks in Jinzhai County and Huangshan District(Table 1)at different segment sizes,including 20,60,and 100 m.The correlation coefficients and RMSE between retrieved ATLAS mean canopy height and corresponding airborne LiDAR metrics were calculated, and used as indicators for assessing the effectiveness of the proposed approach.

3. Results

3.1. Evaluation of the ATL08 products

The comparative analysis of ATL03/ATL08 photons and airborne LiDAR-derived DTM and DSM with different terrain slopes and forest covers (Fig. 5) shows that, overall, the classified photons (ground or canopy)have similar trends,but the canopy photons from ATLAS greatly outnumber the ground photons, indicating the difficulty in extracting DTM from ATLAS data in mountainous regions.Fig.5 also shows that the higher slopes or higher forest covers result in fewer ground photons(Fig.5 a and c).In contrast,relatively low forest covers with gentle slopes generate more ground photons,which are close to airborne LiDAR DTM(Fig.5d),indicating their high data quality.The valley and peak terrain positions considerably influence the quality of ground photons,as shown in Fig.5b and c,where the heights of ground photons are much higher or lower than airborne LiDAR DTM.

Correlation analysis between ATL08 canopy parameters and corresponding parameters from airborne LiDAR(Table 2) showed that,overall,the r values are not high.Most canopy height parameters from ATL08(except for minimum)have higher correlation with the lower percentiles(RH25, RH50, RH60, RH70) and mean from airborne LiDAR, but lower correlation with higher percentiles such as RH80–RHmax. One interesting thing is that Hmin and Hmax from ATL08 have very low r values with corresponding ones from airborne LiDAR, indicating that the maximum or minimum canopy height from the ATLAS ATL08 product are possibly unreliable, at least in this mountainous region. The high uncertainties of Hmin and Hmax may attribute to the outlying noises retained during ATL03 de-noising processing and mislabeled during ATL08 photon classification. Therefore, the scientist team of ATLAS product recommended the 98%canopy height(RH98)be considered as the top of canopy measurement.

Table 2 also shows that mean and median canopy height from ATL08 have the highest r values with those from airborne LiDAR, a conclusion similar to previous studies(e.g.,Chen et al.,2019).In addition,previous studies often used the mean value as the canopy height measurement(Lin et al., 2020; Zhu et al., 2020), thus, it is also used in the followinganalysis.

Table 1 ICESat-2 ATLAS products used in this study.

Fig. 5. A comparison of ATL03/ATL08 photons and airborne LiDAR-derived digital terrain model (DTM) and digital surface model (DSM) with different slopes and forest covers.

Table 2 The correlation between canopy height metrics from ATL08 products and airborne LiDAR based on segment length of 100 m.

The low correlations between ATL08 canopy height parameters and the corresponding ones from airborne LiDAR indicate that the ATL08 canopy heights have high uncertainty. We further examined the correlations between DTM and DSM from ATL08 products and the counterparts from airborne LiDAR and found that although the r values were as high as 0.9997,the RMSE reached 5.9 m for DSM and 8.2 m for DTM.The high uncertainties of both DSM and DTM cause unreliability in the ATL08 product for direct use in biomass estimation or other applications. This situation is especially worse in mountainous regions because of the impacts of steep slopes and dense forest covers on DTM.

3.2. Incorporation of airborne LiDAR DTM and ATLAS data for retrieval of canopy heights

A comparison of ATLAS canopy height calculated based on the combination of airborne LiDAR DTM and ATLAS canopy photon height with the original ATL08 product at segment length of 100 m(see Table 3)indicates considerable improvement reflected by the increase of r and reduction of RMSE when the corresponding airborne canopy height was used as reference.As shown in Table 3,with three strong beams at three acquisition dates, the r values are 0.30–0.56 and RMSE values are 4.49–11.71 m for canopy height of the original ATL08 product.However,they become 0.70–0.94 and 1.32–3.23 m,respectively,after introducing airborne LiDAR DTM, indicating the necessity of using accurate DTM data for retrieval of canopy heights from ATLAS photons.The results also indicate the different degrees of improvement among different ground tracks and acquisition dates, which may be attributed to terrain conditions and forest covers.The improvements are better illustrated in Fig.6,which shows the overestimation of canopy height in the ATL08 product being solved by using airborne LiDAR DTM(Fig.6b).However,there are still some outliers(noises)that need to be removed.This may be achieved by segment filtering proposed in this study.

3.3. Determination of thresholding methods

The performances of different thresholding methods at various segment lengths are presented in Table 4. The result indicates that all three thresholding methods(S, M, and E) increase the r values between ATLAS canopy height and corresponding airborne LiDAR and reduce RMSE compared with non-selection scenario, and the improvement becomes smaller as the segment length becomes longer. Overall, r and RMSE have similar trends:r becomes larger and RMSE becomes smaller as the segment length increases. Although r and RMSE did not have obvious differences among the three thresholding methods,the retention rates have remarkable differences. Therefore, different thresholds determined by thresholding methods should be applied to different segment lengths. Fig. 7 shows great improvement using different thresholding methods compared with no data filtering for different segment lengths.For example,for the short segment lengths of 20 and 40 m,the numbers of total segments are large,and resulting ATLAS canopy heights are highly dispersed(Fig.7a and b);after segment filter using the E method, r increased from 0.69 to 0.75 to 0.79 and 0.82, while RMSE reduced from 2.75 to 3.25 to 2.20 and 2.54 m; although the retention rates are only 21%and 24%,there are still numerous segments retained in the study area.For the segment lengths of 60 and 80 m(Fig.7 c and d),the M method performed better than the S method in removing outliers,while the retention rates were 36.9% and 38.5%, respectively. For the longer segment lengths of over 100 m, the outliers greatly reduced(Fig. 7e–i), and all three thresholding methods had similar r and RMSE values.However,S method kept 76.7%–88.4%of all segments,while M and E methods removed more than half of the segments. Therefore, the optimal thresholding methods for short segment lengths(20 and 40 m),median segment lengths(60 and 80 m),and long segment lengths(over 100 m) are E, M, and S, respectively, which were applied to different ICESat-2 tracks for validation.

3.4. Validation of the proposed canopy height retrieval and filtering method

According to the thresholding methods selected in section 3.3 for different segment lengths,i.e.,E for 20 and 40 m,M for 60 and 80 m,and S for over 100 m, we identified the optimal thresholds, calculated the mean canopy heights through a combination of ATLAS photon attributes and airborne LiDAR DTM at segment lengths of 20, 60, and 100 m for different ground tracks, and compared them with the corresponding airborne LiDAR canopy heights in terms of r and RMSE (Table 5). We found that for most tracks from all acquisition dates, except for GT3 of Jun 10, 2019 and GT1 of August 20, 2019, r values between retrieved ATLAS canopy height and reference height from airborne LiDAR achieve 0.70, and RMSE are lower than 3.00 m, indicating the retrieved ATLAS canopy heights at different segment sizes are much close to the reference canopy height from airborne LiDAR, implying the proposed method is effective for canopy height retrieval from ATLAS data.We also found that for the same track,increase of segment length improves the result quality in terms of increased r and reduced RMSE values. For different ground tracks of the same date, the accuracies vary; for instance, r for GT1 on August 20,2019 at segment length of 20 m is only 0.61 and RMSE is 4.35 m,but for GT2 and GT3 on the same date,r values are 0.89 and 0.91 and RMSE values are 2.05 and 1.90 m, respectively, although the threshold for GT1 is greater than the other two. Even using the same threshold, r and RMSE values are much different among tracks, exemplified by GT1 and GT3 at a 20-m segment on September 9,2019 and 60-m segment on August 20,2019.This also indicates that when r reaches similar values,the required thresholds vary, depending on the data acquisition conditions and geophysical characteristics,such as topography and vegetation cover.

From Table 5 we also found that after segment filtering using the thresholds determined by thresholding methods optimal for different segment lengths, the retention rate for 20, 60, and 100 m were 20.1%–23.8%, 34.4%–59.5%, and 74.4%–85.6%, respectively. Therefore, we took the lowest retention rates for corresponding segment lengths,that is 20%,35%,and 75%as an alternative approach to determine thresholds for segments 20, 60, and 100 m. This would ensure high quality of retained segments. The accuracy assessment results are presented in Table 6. Comparative analysis of the contents of Tables 5 and 6 shows similar r and RMSE results, indicating that we can directly use the proposed retention rate to retrieve canopy heights, which is much simpler and easier.

4. Discussion

4.1. Impacts of topography and forest covers on quality of the retrieved canopy heights

Previous studies using ICESat-2 ATLAS data are mainly located ingentle terrain (Li et al., 2020). Rarely has research explored the mountainous regions with dense forest covers, because very low numbers of ground photons do not generate good-quality DTM data. Previous research indicated that topography and vegetation cover are important factors influencing the accuracy of terrain elevation retrieval (Wang et al., 2019; Neuenschwander et al., 2020). Our research in subtropical mountainous regions indicates the RMSE of 8.09 m for DTM at segment length of 100 m between the ATLAS product and airborne LiDAR DTM.This is because the ATLAS products used global multi-resolution terrain elevation data as DTM, which has relatively high errors; for example,RMSE can be 26–30 m at segment length of 250 m(Danielson and Gesch,2011). Another problem is that the ATLAS product cannot effectively extract DTM at valley and peak points (Neuenschwander et al., 2019b),resulting in high uncertainty in undulating terrain,in turn leading to high uncertainty in canopy height retrieval. Our evaluation of the ATL08 product with airborne LiDAR in subtropical mountainous regions indicated when slopes are ＜15°, 15°–30°, and 30°–45°, the correlation coefficients between mean canopy heights of ATL08 and airborne LiDAR reduced from 0.60, to 0.47, then to 0.10, respectively, while the RMSE increased from 7.02 to 8.56 m,then to 9.8 m,respectively,implying the accuracy of canopy height decreases as slope increases. Therefore, this research proposed a combination of high-resolution airborne LiDAR DTM and ATLAS photon heights to generate canopy height at various segment lengths, which provided much improved accuracy of canopy height products. Because DTM usually does not change over time and precise DTM from airborne LiDAR are available at large scale in many places,such as Fujian Province in 2013–2016 and Anhui Province in 2015–2016, incorporation of this kind of DTM data with 2019–2021 ATLAS products can quickly and accurately update canopy height products,providing new data sources for forest stock volume or biomass modeling in a large area.

Table 3 A comparison of canopy height assessment results with and without use of the airborne LiDAR digital terrain model (DTM) based on segment length of 100 m.

Fig. 6. A comparison of mean canopy height with and without use of airborne LiDAR digital terrain model(DTM)based on the ground track GT1L data on June 10,2019(r,correlation coefficient;RMSE,root mean square error).(a)Canopy height from ATL08,(b)Canopy height based on combination of ATLAS photon height and airborne LiDAR DTM.

Table 4 A comparison of correlation coefficient(r),root mean square error(RMSE)(m),and retention rate(%)based on different segment lengths,thresholding methods,and thresholds.

4.2. Determination of suitable methods for selection of reliable segments

Fig. 7. Comparison of different thresholding methods at various segment lengths. (a)–(i) represent the segment lengths from 20 to 180 m at intervals of 20 m; the color dots represent retention segments after using three thresholding methods; r and RMSE are based on the selected thresholding methods for associated segment lengths:(a)and(b)from E,(c)and(d)from M,and(e)–(i)from S.(For interpretation of the references to color in this figure legend,the reader is referred to the Web version of this article.)

Removal of unreliable segments is needed for retrieval of accurate canopy height from ATLAS data.Previous research has explored different approaches to filter the unreliable segments: (1) removal of 30-m grids when the number of ground photons is less than 4 in a grid cell(Lin et al.,2020) to ensure high accuracy of DTM; (2) removal of 30-m grid cells when the number of photons within the unit is less than 40 or greater than 170,securing the accuracy of extracted canopy heights(Zhu et al.,2020);(3)comprehensively considering the impacts of the data retention rate and different thresholds on the retrieval of forest canopy height,for example, keeping grids that contain at least 3 segments of an ATL08 product within 1-km grid(Sun et al.,2020).In reality,the mountainous regions with steep slopes and dense forest covers have much fewer ground photons,as shown in Fig.1,barely producing DTM data.Because very little research had explored the filtering methods for removing unreliable segments in mountainous regions,this research explored three thresholding methods used to determine the thresholds for different segment lengths,based on which lower-quality segments were removed.However, the number of photons that return to the detection devices depends on surface reflectivity and atmospheric conditions, and the optimal threshold of photons at a specific segment length identified from the three methods vary among different ground tracks and different acquisition dates. Generally, the E thresholding method is suitable for finer scale (20 and 40 m), the M thresholding method is suitable for medium scale (segment length 40–60 m), while the S thresholding method is suitable for coarser scale(segment length ≥100 m).With these methods and segment lengths, both the accuracies of retrieved canopy height of retained segments and the numbers of retained data are satisfied. From the experiment conducted in the study area, we proposed a simple approach based on the retention rate for different segment lengths and obtained reliable canopy height data,and our approach was further proven transferable after examining different tracks. Specifically, the retention rates are 20%, 35%, and 75% for segment lengths of 20–40,60–80 m,and greater than 100 m, respectively.

4.3. Impacts of segment lengths on the retrieval of canopy height

Although some previous research explored the scale issue related to canopy height retrieval (Chen et al., 2019; Huang et al., 2020; Li et al.,2020), no consistent conclusions were obtained; in particular, noresearch was done in subtropical mountainous regions with complex steep slopes and dense forest covers. A simulation study using airborne LiDAR data to explore retrieval of canopy height from ICESat-2 at different spatial resolutions (10, 16, 20, 30, 40, and 50 m) found R2(coefficient of determination)of 0.69 and RMSE of 2.2 m at 50 m(Chen et al., 2019). Huang (2021) explored the ATLAS data for retrieving canopy parameters with different window sizes from 20 to 50 m at 5-m intervals, and found the best window size of 20 m with RMSE of 4.28 m.Another study in Inner Mongolia for canopy height retrieval at spatial resolutions of 10,30,250,500,and 1000 m found that the mean canopy height at 250 m performed best with RMSE of 2.39 m (Li et al., 2020).These studies conducted in regions with relatively gentle slopes did not provide general conclusions.No research has explored the spatial issues for canopy height retrieval in mountainous regions with complex topographic conditions and dense forest covers due to the lack of precise DTM data.

Table 5 Comparison of accuracy of the retrieved canopy heights based on selected thresholding methods (E for 20 m, M for 60 m, and S for 100 m).

Table 6 Comparison of accuracy assessment results of the retrieved canopy heights based on determined retention rates(20%for 20 m,35%for 60 m,and 75%for 100 m).

Our research explored different scales from 20 to 200 m at 20-m intervals and provides filtering approaches corresponding to the scales.The considerably improved r values and reduced RMSE values (see Table 6)indicate the effectiveness of the proposed method.The research also provides guidelines to select suitable scales for linking the canopy parameters to different spatial resolution images such as Sentinel-2,Landsat, and MODIS data for wall-to-wall mapping of forest canopy height.

5. Conclusions

The ICESat-2 ATLAS provides new opportunity for retrieval of canopy heights at national and global scales,but the terrain conditions and forest covers in mountainous regions result in high uncertainty in canopy height retrieval. We propose incorporation of airborne LiDAR DTM and ATLAS product to generate canopy heights at different segment lengths,which provides much improved results.Specifically,the ATLAS product in this study area has relatively poor correlation with airborne LiDARderived canopy height; that is, the r values are between 0.30 and 0.53,and RMSE are between 4.49 and 11.71 m, whereas by incorporating airborne LiDAR DTM and ATLAS data, the r values increased to 0.70–0.94, and RMSE decreased to only 1.32–3.23 m. This research found that the retention rate can be used to select high-quality segments at various lengths in mountainous regions:the retention rate of 20%for segment lengths of 20 and 40 m;35%for 60 and 80 m,and 75%for 100 m and over. This method can ensure the selected segments effectively represent DSM, while the combination of airborne LiDAR DTM and ATLAS products can ensure the retrieved canopy heights have high r and low RMSE values compared to LiDAR canopy heights. According to evaluation of canopy height results based on 20,60,and 100 m segments,the r values of 0.61–0.91,0.65–0.92,and 0.68–0.94,respectively,can be reached, and RMSE can be as low as 1.90–4.35, 1.55–3.63, and 1.34–3.23 m, respectively. Because high precision DTM data are available at a large scale, the proposed method provides an easy way to quickly retrieve canopy heights from ICESat-2 ATLAS products in mountainous regions with complex terrain conditions and dense forest covers.

Funding

This study was financially supported by the National Natural Science Foundation of China(No.32171787).

Declaration of competing interest

The authors declare no conflicts of interest.