Correction of Terrain Effects on Forest Canopy Height Estimation Using ICESat-2 and High Spatial Resolution Images

: The Ice, Cloud, and Land Elevation Satellite–2 (ICESat–2) carries the Advanced Topographic Laser Altimeter System (ATLAS), enabling global canopy height measurements from forest canopy height models (CHMs). Topographic slope is a crucial factor affecting the accuracy of canopy height estimates from ICESat–2 CHMs, but it has not been sufﬁciently studied. This paper aims to eliminate the inﬂuence of slope on canopy height estimates from ICESat–2 data and establishes a method for correcting forest canopy heights based on high spatial resolution digital orthophoto maps (DOM). The cross-track photons are corrected horizontally to eliminate the estimation error. Multi-resolution segmentation is used to segment tree crowns in the DOM, and the distance and relative position between the top of canopy (TOC) photons and the center point of the crown are calculated. TOC photon correction rules are established for different terrains, and the vertical error of the TOC photons is corrected. The results indicate that the vertical error increases exponentially with the slope. The cross-track photon correction and the TOC photon correction method eliminate the effect of slope on canopy height estimates. The cross-track photon correction method reduces the mean absolute error (MAE) and root mean square error (RMSE) of the canopy height estimates by 35.71% and 35.98%, respectively. The TOC photon correction approach further reduces the MAE and RMSE by 23% and 19.23%, respectively. The proposed method has signiﬁcantly higher accuracy for forest canopy height estimation using ICESat–2 data than the traditional method.


Introduction
Vegetation canopy height is a critical input parameter for vegetation biomass modeling, which is crucial for tracking the spatial and temporal dynamics of forest ecosystems [1][2][3] and providing information on how the carbon cycle and terrestrial ecosystems affect the future climate [4][5][6]. Airborne laser scanning (ALS) is a relatively mature method widely used in modern forest monitoring and has replaced fieldwork in many cases [7][8][9]. Although a growing number of countries are acquiring ALS data, collecting these data over large areas is extremely expensive. Current opportunities for obtaining large-scale LiDAR data are provided by the Global Ecosystem Dynamics Investigation (GEDI) project and the Ice, Cloud, and Land Elevation Satellite-2 (ICESat-2) [10][11][12].
The Advanced Topographic Laser Altimeter System (ATLAS) onboard ICESat-2 is a single-photon sensitive LiDAR system capable of detecting individual photons from reflected pulses [13][14][15]. ICESat-2 has 4 levels of 22 products, ATL00~ATL21. ATL03 and ATL08 are digital elevation model (DTM) and digital surface model (DSM) products. ATL03 is the global geolocated photon data, and ATL08 is the land and vegetation height data. The latter is a collection of ground and canopy elevations obtained by filtering the the authors did not report the source of the slope error and whether it could be eliminated. Determining the errors in forest height estimates due to slope in complex topographic areas has become a critical issue for using ICESat-2 data for large-scale biomass estimation. A methodology to attenuate the influence of surface terrain on canopy height estimates is urgently needed.
The objective of this study is to propose a physical geometric correction model based on high spatial resolution digital orthophoto map (DOM) data to correct canopy height estimates and improve the accuracy of ICESat-2 canopy parameter extraction in complex terrain. We focus on the following aspects: (1) we establish an ICESat-2 canopy height correction method that eliminates the effect of slope to enable ICESat-2 canopy height inversion. (2) We investigate the accuracy of ATL03 and ATL08 canopy data in subtropical regions and the effect of different window sizes and percentages of canopy height parameters on forest canopy height estimation. (3) The effects of the photon offset distance, canopy relief rate, and canopy density on the inversion accuracy of the canopy height estimates derived from the point cloud are analyzed. The results of this study provide new methods and ideas for ICESat-2 canopy data extraction in complex terrain.

Study Area
The study area is a forested region in Nanning City, Guangxi Zhuang Autonomous Region, China (22.57 • N~23.00 • N, 108.23 • E~108.24 • E (cf. Figure 1). Nanning City has a subtropical monsoon climate, with an annual average temperature of about 21.6 • C and average annual precipitation of 1241~1753 mm [36]. The study area is hilly land with an elevation range of 120~360 m (World Geodetic System 1984, WGS84) and a slope range of 0 •~5 1 • , with slope greater than 15 • accounting for 64.46% of the study area. The study area has coniferous, broadleaf, and mixed coniferous forests. The forest cover exceeds 92%, and the canopy height has a range of 0.1~24.6 m. The average tree density is 1300 trees per hectare, the average stand basal area is 11.30 m 2 /hm 2 , and the average stand volume is 147.47 m 3 .

ICESat-2 Data
ICESat-2 has an average orbital altitude of 496 km (±6 km) and covers an area ranging from 88°N to 88°S [30,37,38]. The footprint of ATLAS has a diameter of about 14 m and a footprint spacing of 0.7 m along the track [39]. ATLAS uses three pairs of beams, each pair consisting of a laser beam with intensity-to-weak beam energy ratio of 4:1 [12] (cf. Figure  2). During the light pulse emission, the photon triggers a timer as it passes through the beam splitter and when it returns and passes through the filter. The time-stamped data of the returning photons are transmitted to the ground via the electronics and communication system on ICESat-2. The WGS84 ellipsoidal heights of the surface are determined by calculating the position of the spacecraft and the distance traveled by the photon (https://icesat-2.gsfc.nasa.gov/space-lasers, accessed on 20 October 2021). This study only analyzed the canopy height of strong nighttime beams. The track ID numbers are ATL03_20191003151649_01030501_003_01_gt3r and ATL08_2019100 3151649_01030501 _003_01_gt3r, and the data were obtained on 3 October 2019 from https://search.earthdata. nasa.gov (accessed on 5 March 2021).

Research Data 2.2.1. ICESat-2 Data
ICESat-2 has an average orbital altitude of 496 km (±6 km) and covers an area ranging from 88 • N to 88 • S [30,37,38]. The footprint of ATLAS has a diameter of about 14 m and a footprint spacing of 0.7 m along the track [39]. ATLAS uses three pairs of beams, each pair consisting of a laser beam with intensity-to-weak beam energy ratio of 4:1 [12] (cf. Figure 2). During the light pulse emission, the photon triggers a timer as it passes through the beam splitter and when it returns and passes through the filter. The time-

ALOS PALSAR DEM Data
The Advanced Land Observing Satellite (ALOS) has a Phased Array type L-band Synthetic Aperture Radar (PALSAR) capable of almost all-weather observations [40]. The ALOS PALSAR DEM has a pixel size of 12.5 m [41]. The file (AP\U 16054\U FBS\U F0440\U RT1) was downloaded from https://asf.alaska.edu (accessed on 20 March 2021) and contained data obtained in October 2018. Due to the relatively high pixel resolution and accuracy of ALOS PALSAR DEM data, it is used as the reference elevation data.

CAF-LiCHy LiDAR Data
A LiDAR sensor with a small spot size is preferable for evaluating the spatial structure of forests [42][43][44]. The Institute of Resource Information at the Chinese Academy of Forestry has developed an integrated remote sensing observation platform called CAF-LiCHy (Chinese Academy of Forestry, LiDAR, CCD, Hyperspectral). The system consists of a full-waveform airborne lidar (Riegl LMS-Q680i) system, an airborne push-broom hyperspectral scanner (AISA Eagle II), and a high-resolution charge-coupled device (CCD) camera (DigiCAM-H60) [15,45]. The CAF-LiCHy is similar to the LiDAR Hyperspectral Thermal (G-LiHT) system of the National Aeronautics and Space Administration (NASA) Goddard Space Flight Center. The CHM extracted from the CAF-LiCHy LiDAR data was used to verify the accuracy of the canopy height model (CHM) obtained from the ATLAS data. The spatial resolution and vertical accuracy of the CAF-LiCHy data are 1.0 m and 0.15 m, respectively (cf. Table 1). We obtained data in April 2018 in the WGS-84 coordinate system.

ALOS PALSAR DEM Data
The Advanced Land Observing Satellite (ALOS) has a Phased Array type L-band Synthetic Aperture Radar (PALSAR) capable of almost all-weather observations [40]. The ALOS PALSAR DEM has a pixel size of 12.5 m [41]. The file (AP\U 16054\U FBS\U F0440\U RT1) was downloaded from https://asf.alaska.edu (accessed on 20 March 2021) and contained data obtained in October 2018. Due to the relatively high pixel resolution and accuracy of ALOS PALSAR DEM data, it is used as the reference elevation data.

CAF-LiCHy LiDAR Data
A LiDAR sensor with a small spot size is preferable for evaluating the spatial structure of forests [42][43][44]. The Institute of Resource Information at the Chinese Academy of Forestry has developed an integrated remote sensing observation platform called CAF-LiCHy (Chinese Academy of Forestry, LiDAR, CCD, Hyperspectral). The system consists of a full-waveform airborne lidar (Riegl LMS-Q680i) system, an airborne push-broom hyperspectral scanner (AISA Eagle II), and a high-resolution charge-coupled device (CCD) camera (DigiCAM-H60) [15,45]. The CAF-LiCHy is similar to the LiDAR Hyperspectral Thermal (G-LiHT) system of the National Aeronautics and Space Administration (NASA) Goddard Space Flight Center. The CHM extracted from the CAF-LiCHy LiDAR data was used to verify the accuracy of the canopy height model (CHM) obtained from the ATLAS data. The spatial resolution and vertical accuracy of the CAF-LiCHy data are 1.0 m and 0.15 m, respectively (cf. Table 1). We obtained data in April 2018 in the WGS-84 coordinate system. Ancillary image data included DOM data from the CAF-LiCHy system and highresolution satellite image data from Google Earth. The DOM data were acquired by the CCD digital camera (DigiCAM-H60), which has a focal length of 50 mm. The images have 60 million pixels and 0.2 m spatial resolution [45]. The data were obtained in April 2018. Table 1 list the sensors' parameters. The satellite image data from Google Earth were L19 level (spatial resolution of approx. 0.3 m) and were acquired in March 2018 and accessed on 7 August 2022. The image data had a GeoTIFF file format and were in the WGS-84 vertical coordinate system.

Methodologies
The estimation of canopy height parameters requires explicit TOC and ground photon signals. In this study, the ICESat-2 data processing included cross-track photon correction, photon denoising, photon classification, TOC photon correction, canopy parameter extraction, and accuracy evaluation. The canopy heights obtained from the ATL08 data (CHM ATL08 ), ATL03 original data (CHM ATL03 initial ), ATL03 data after cross-track photon correction (CHM CTPC ), ATL03 data after cross-track photon correction and Google Image data after canopy height correction (CHM CCR-1 ), and ATL03 data after cross-track photon correction and airborne image canopy height correction (CHM CCR-2 ) were compared with the CHM ALS extracted from the CAF-LiCHy LiDAR Data. Figure 3 depicts the flowchart of the inversion of different CHMs. Ancillary image data included DOM data from the CAF-LiCHy system and highresolution satellite image data from Google Earth. The DOM data were acquired by the CCD digital camera (DigiCAM-H60), which has a focal length of 50 mm. The images have 60 million pixels and 0.2 m spatial resolution [45]. The data were obtained in April 2018. Table 1 list the sensors' parameters. The satellite image data from Google Earth were L19 level (spatial resolution of approx. 0.3 m) and were acquired in March 2018 and accessed on 7 August 2022. The image data had a GeoTIFF file format and were in the WGS-84 vertical coordinate system.

Methodologies
The estimation of canopy height parameters requires explicit TOC and ground photon signals. In this study, the ICESat-2 data processing included cross-track photon correction, photon denoising, photon classification, TOC photon correction, canopy parameter extraction, and accuracy evaluation. The canopy heights obtained from the ATL08 data (CHM ATL08), ATL03 original data (CHM ATL03 initial), ATL03 data after crosstrack photon correction (CHM CTPC), ATL03 data after cross-track photon correction and Google Image data after canopy height correction (CHM CCR-1), and ATL03 data after cross-track photon correction and airborne image canopy height correction (CHM  were compared with the CHMALS extracted from the CAF-LiCHy LiDAR Data. Figure 3 depicts the flowchart of the inversion of different CHMs.    [40,46], and the distance of the photon's cross-track on flat land or on low slopes does not affect the ground elevation and canopy height estimation. However, on large slopes and complex terrain, disregarding the cross-track offset of the photons might cause errors in the DTM and DSM (cf. Figure 4a). Figure 4d shows an irregular distribution of photons on each side of the track, indicating that the cross-track offsets are not the same. Therefore, before the photon data were transformed into a photon point cloud, photon correction was performed using the cross-track photon correction method proposed by Li et al. [47]. First, we obtained the mapping points of the cross-track photons (Equation (1)). The reference elevations of the original photon (P) and the mapping point (P') of ATL03 were extracted from the ALOS PALSAR DEM (a cubic interpolation method was used to interpolate the DEM to 1.0 m spatial resolution) as the reference elevations z 0 and z 1 of the original photon and the mapping point, respectively. Photons with a cross-track distance ≥0.50 m or ≤−0.50 m were corrected to the reference ground track (cf. Figure 4c), and the corrected photon elevation was obtained from Equation (2).
where x 0 , y 0 are the original latitude and longitude of the ATL03 photon, respectively; x 1 , y 1 are the latitude and longitude of the mapping point, respectively; z 0 , z 1 are the original elevations of the photon and mapping point, respectively, E C is the elevation of the mapping point, and E i is the original elevation of the ATL03 photon.  [40,46], and the distance of the photon's cross-track on flat land or on low slopes does not affect the ground elevation and canopy height estimation. However, on large slopes and complex terrain, disregarding the cross-track offset of the photons might cause errors in the DTM and DSM (cf. Figure 4a). Figure 4d shows an irregular distribution of photons on each side of the track, indicating that the cross-track offsets are not the same. Therefore, before the photon data were transformed into a photon point cloud, photon correction was performed using the cross-track photon correction method proposed by Li et al. [47]. First, we obtained the mapping points of the cross-track photons (Equation (1)). The reference elevations of the original photon (P) and the mapping point (P') of ATL03 were extracted from the ALOS PALSAR DEM (a cubic interpolation method was used to interpolate the DEM to 1.0 m spatial resolution) as the reference elevations z0 and z1 of the original photon and the mapping point, respectively. Photons with a cross-track distance ≥0.50 m or ≤−0.50 m were corrected to the reference ground track (cf. Figure 4c), and the corrected photon elevation was obtained from Equation (2).
(1) (2) where x0, y0 are the original latitude and longitude of the ATL03 photon, respectively; x1, y1 are the latitude and longitude of the mapping point, respectively; z0, z1 are the original elevations of the photon and mapping point, respectively, EC is the elevation of the mapping point, and Ei is the original elevation of the ATL03 photon.

Noise Photon Removal
The photon's geographic location was converted to the along-track distance (ATD) before noise removal. The data were first coarsely denoised to narrow down the range of possible signal photons. The photon cloud data were divided into different along-track photon cloud data using a window size of 20 m in the ATD. Histograms of the elevation statistics were calculated, and the histogram peak positions were obtained for each window. They were used as the center position of the signal photons, and an elevation threshold was set near the center position to establish a buffer of photon signals (Equations (3) and (4)). Subsequently, we removed the noise photons outside the buffer.
where minr is the minimum elevation, maxr is the maximum elevation, h α is the elevation at the histogram peak position, and σ is the standard deviation of the photon elevation within each window. The improved Ordering Points to Identify the Clustering Structure (OPTICS) algorithm proposed by Zhu et al. [17] was used for fine denoising. The improved OPTICS algorithm consists of three steps. First, an elliptical search region was established. Second, the core distance and reachable distance of each photon were calculated based on the neighborhood radius and neighborhood minimum number of points. Third, the photons were sorted according to the reachable distance from smallest to largest, and the optimal reachable distance threshold of each window was determined using the maximum inter-class variance algorithm. Photons larger than the threshold were regarded as noise points, and the others were regarded as signal points. Figure 5a shows the ATL03 cross-track correction data and the result after noise removal.

Photon Classification
Following photon noise removal, the data contained few ground photons and numerous TOC photons, making it difficult to distinguish them and accurately depict the terrain and canopy [48]. An overlapping moving window was created to identify the ground and canopy photons.
We provide a brief summary of the methods used to classify ground photons. They include overlapping moving windows, the percentile method, and the slope filtering method. First, a window was created. A correlation exists between the extraction error and the extraction window size and step size [40,41,49]. Preliminary investigations indicated that the moving window width to detect ground photons should be 20 m, and the moving step length should be 10 m to enable the precise extraction of the photons of interest. Second, an elevation percentile range was selected to extract potential ground photons within the window range. Since the ATL03 data in this region had a low photon density, photons in the elevation percentile range of [0.00, 0.15] were chosen to ensure that a sufficient number of photons were available for terrain height estimation. However, near-ground noise photons and near-ground vegetation photons remained in the photon point cloud after denoising, and the photons in the window elevation percentile range of [0.00, 0.15] were not necessarily all ground photons. Ground photons have a continuous distribution, whereas near-ground noise photons and near-ground vegetation photons lack continuity compared with the elevation values of neighboring points. These photons are prone to sudden elevation changes (sudden increase or decrease) because the data have not been smoothed. The slope filtering method (SFMB) proposed by Li et al. [47] was used to minimize the effect of near-ground photon noise on ground photon extraction and extract possible ground photons obtained in the previous step. The SFMB method consists of two steps.  The slope of the possible ground photons was calculated, and the maximum slope in the reference elevation data ALOS PALSAR DEM was used as the threshold value. Photons larger than the threshold value were marked as noise photons. The maximum value of the slope change in the reference elevation data was used as the threshold value, and photons smaller than the threshold value were marked as ground photons. Figure 5b depicts the results (red dots).
The correct extraction of TOC photons may be hampered by photon noise in the photon cloud remaining after denoising [50]. We used the TOC classification method proposed by Popescu et al. [20] to determine the TOC photons. The width of the moving window and the movement steps were the same as for the ground photon classification. The TOC noise photons were those in the given elevation percentile range [0.99, 1] for each moving window. After removing the photons in this percentile range, the elevation percentile was recalculated, and the photon cloud corresponding to the photon elevation percentile range of [0.96, 0.99] was extracted and used as TOC photons. The canopy photon results are shown as green dots in Figure 5c.

TOC Photon Correction Effect of Slope
The canopy height photon point cloud contains information on the height difference between the DSM and DTM on the same coordinate plane. Although a photon-counting LIDAR system has the advantages of providing high frequency and high-resolution data, the TOC photons are not necessarily located at the top of the canopy, resulting in errors in the canopy height estimation. The difference between the elevation at the TOC and that on the ground should be the canopy height of trees. However, this height difference differs from the real canopy height on slopes (such as on hills, valleys, steep slopes, or degraded places) [51].
When the difference between canopy photon data and DTM was used to compute canopy height, Figure 6 illustrates the impact of the surface slope on canopy height based on the difference between the canopy height of the photon data and the DTM. The LiDAR variable RH100 is defined as the distance between the apex of the canopy and the ground point. On flat land, RH100 is approximately the vegetation canopy height (cf. Figure 6a). When a tree is located on a slope, the photons hit the same position in the canopy, and the difference between the canopy photons and the DTM is h a '. This parameter has been used as the canopy height in most studies on canopy height inversion based on ICESat-2 data. However, the actual canopy height is the difference between the height of the TOC photons and the height at the point where the tree roots enter the ground, i.e., h a in Figure 6b. The difference (or error) between h a and h a ' is caused by the slope and is defined by Equation (5). The position of the photon relative to the trunk center point affects how the slope error is corrected. Figure 6b,c show this difference at different slope heights. The height discrepancy between causes errors in estimating the canopy height and is largely terrain-induced. Therefore, the effect of topography must be eliminated to extract the canopy height accurately.
The distance between the photon location and the trunk center point (D d ), the relative position, and the slope affect the terrain-induced canopy height error. Table 2 lists the maximum errors for different slopes and different D d . The vertical displacement caused by the slope increases exponentially with the slope of the terrain. The terrain has a significant impact on canopy height measurements that cannot be ignored in extremely steep terrain. The distance between the TOC photon location and the center point of the crown cannot be directly obtained from probability sampling of ATLAS photons [24,52,53]. Therefore, we utilize crown segmentation to quantify the effect of the distance between the TOC photon location and the center point of the crown on canopy height extraction.
where h is the tree height obtained by inversion, h id actual tree height, D d is the horizontal distance between the TOC photon location and the center point of the crown, and p is the slope.
difference between the canopy photons and the DTM is ha'. This parameter has been used as the canopy height in most studies on canopy height inversion based on ICESat-2 data. However, the actual canopy height is the difference between the height of the TOC photons and the height at the point where the tree roots enter the ground, i.e., ha in Figure  6b. The difference (or error) between ha and ha' is caused by the slope and is defined by Equation (5). The position of the photon relative to the trunk center point affects how the slope error is corrected. Figure 6b,c show this difference at different slope heights. The height discrepancy between causes errors in estimating the canopy height and is largely terrain-induced. Therefore, the effect of topography must be eliminated to extract the canopy height accurately. The distance between the photon location and the trunk center point (Dd), the relative position, and the slope affect the terrain-induced canopy height error. Table 2 lists the maximum errors for different slopes and different Dd. The vertical displacement caused by the slope increases exponentially with the slope of the terrain. The terrain has a significant impact on canopy height measurements that cannot be ignored in extremely steep terrain. The distance between the TOC photon location and the center point of the crown cannot be directly obtained from probability sampling of ATLAS photons [24,52,53]. Therefore, we utilize crown segmentation to quantify the effect of the distance between the TOC photon location and the center point of the crown on canopy height extraction.
where ℎ is the tree height obtained by inversion, ℎ id actual tree height, is the horizontal distance between the TOC photon location and the center point of the crown, and is the slope.  Crown Segmentation Multi-resolution segmentation (MS) is commonly used for segmenting ultra-high resolution images [54][55][56]. We used eCognition software (Definens Developer 8.7) for crown segmentation of the DOM with a resolution of 0.2 m and Google Earth images with a resolution of 0.3 m. A 9 × 9 median filter was applied to the DOM prior to segmentation to prevent over-segmentation [51]. The scale, shape ratio, and compactness ratio are required input parameters for MS in eCognition. The scale parameter affects the average size of the segmented objects, the shape ratio controls the homogeneity, and the compactness ratio influences the object's shape. The optimal segmentation parameters for various forest types were determined after extensive experiments and considering the vegetation structure and distribution. The 0.2 m resolution image had a scale parameter of 12~16, a shape ratio of 0.7~0.8, and a compactness ratio of 0.6~0.9. The 0.3 m resolution image had a scale parameter range of 10~12, a shape ratio range of 0.4~0.5, and a compactness ratio range of 0.6~0.7. The To assure quality, the final segmentation results were visually personally examined for an appropriate segmentation performance. The segmentation results are shown in Figure 7.

TOC Photon Correction
The slope and slope direction parameters were calculated from the DTM derived from the ICESat-2 data. The relative positions of the TOC photons and crown center points were obtained after segmenting the high spatial resolution DOM. We used the slope, Dd, the orientation parameter, and the canopy correction rule (CCR) to eliminate the canopy height errors caused by the slope (cf. Equation (7)). If the TOC photon and the crown center point were located on the ATD, the relative distance could be calculated directly. Otherwise, the crown center point was projected onto the ATD, and the relative distance was calculated. The TOC photon locations whose Dd were greater than 0.5 m were corrected. The values were not corrected when the TOC photons were located between adjacent canopies (cf. Figure 8). The highest local canopy point was obtained by Equation (7) and represented the corrected canopy height. Due to the different crown segmentation results for the 0.2 m and 0.3 m spatial resolution images, the locations of the crown centroids were different, resulting in differences in Dd and S. The correction results obtained from the 0.3 m spatial resolution image were called CCR-1, and those based on the 0.2 m spatial resolution image were called CCR-2.

ℎ
• tan where ℎ is the canopy height after correction, is the TOC photon height, Dd is the distance between the TOC photon and the center point of the crown, is the slope, and S is the orientation parameter. When the crown center point was located on the south slope, the TOC photon was located south relative to the crown center at S = 1. The TOC photon was located north relative to the crown center at S = −1. When the crown center point was located on the north slope, the TOC photon was located south relative to the crown center at S = −1, and the TOC photon was located north relative to the crown center at S = 1. When Dd ≤ 0.5 m or the TOC photon was not located on the tree crown, S = 0. The calculations were performed in ArcGIS 10.7 and Python 3.8. The geometric moment of the images was calculated to derive a binary map and obtain the center of mass position of the segmented object (Equation (6)), which was used as the crown center [57]. We vectorized the final segmentation results to obtain the tree crown delineations.
where φN(x, y) is the pixel value of the binary image, C 10 is the cumulative sum of the product of the x-coordinate value of the tree canopy profile and the gray value of the image element, C 01 is the cumulative sum of the product of the y-coordinate value of the tree canopy profile and the gray value of the image element, and C 00 is the cumulative sum of the gray values of the canopy profile image elements.

TOC Photon Correction
The slope and slope direction parameters were calculated from the DTM derived from the ICESat-2 data. The relative positions of the TOC photons and crown center points were obtained after segmenting the high spatial resolution DOM. We used the slope, D d , the orientation parameter, and the canopy correction rule (CCR) to eliminate the canopy height errors caused by the slope (cf. Equation (7)). If the TOC photon and the crown center point were located on the ATD, the relative distance could be calculated directly. Otherwise, the crown center point was projected onto the ATD, and the relative distance was calculated. The TOC photon locations whose D d were greater than 0.5 m were corrected. The values were not corrected when the TOC photons were located between adjacent canopies (cf. Figure 8). The highest local canopy point was obtained by Equation (7) and represented the corrected canopy height. Due to the different crown segmentation results for the 0.2 m and 0.3 m spatial resolution images, the locations of the crown centroids were different, resulting in differences in D d and S. The correction results obtained from the 0.3 m spatial resolution image were called CCR-1, and those based on the 0.2 m spatial resolution image were called CCR-2.
h ATLCCM = H TOC − S·D d tan p where h ATLCCM is the canopy height after correction, H TOC is the TOC photon height, D d is the distance between the TOC photon and the center point of the crown, p is the slope, and S is the orientation parameter. When the crown center point was located on the south slope, the TOC photon was located south relative to the crown center at S = 1. The TOC photon was located north relative to the crown center at S = −1. When the crown center point was located on the north slope, the TOC photon was located south relative to the crown center at S = −1, and the TOC photon was located north relative to the crown center at S = 1. When D d ≤ 0.5 m or the TOC photon was not located on the tree crown, S = 0. The calculations were performed in ArcGIS 10.7 and Python 3.8.

Extraction of Canopy Parameters
The DSM and DTM were obtained by applying a third spline interpolation to the ground photons and TOC photons (cf. Figure 5d). The DTM and DSM were subtracted to obtain the CHM. The CHM was divided into 10 m, 20 m, 30 m, 40 m, 50 m, 60 m, 70 m, 80 m, 90 m, and 100 m segments in the along-track direction, and the effect of the segment size on the accuracy was compared. We calculated the relative height index for the different segments (RH70, RH75, RH80, RH85, RH90, RH95, RH98, RH100). In addition, we used the canopy relief ratio (CRR) [58,59] to analyze the canopy height differences within a 10 m segment: where is the height of each vegetation point relative to the ground, is the canopy elevation, is the ground elevation, DH is the dataset obtained by sorting the height of all vegetation points relative to the ground in ascending order, i is the

Extraction of Canopy Parameters
The DSM and DTM were obtained by applying a third spline interpolation to the ground photons and TOC photons (cf. Figure 5d). The DTM and DSM were subtracted to obtain the CHM. The CHM was divided into 10 m, 20 m, 30 m, 40 m, 50 m, 60 m, 70 m, 80 m, 90 m, and 100 m segments in the along-track direction, and the effect of the segment size on the accuracy was compared. We calculated the relative height index for the different segments (RH70, RH75, RH80, RH85, RH90, RH95, RH98, RH100). In addition, we used the canopy relief ratio (CRR) [58,59] to analyze the canopy height differences within a 10 m segment:  (11) where HC i is the height of each vegetation point relative to the ground, H canopy is the canopy elevation, H terrain is the ground elevation, DH is the dataset obtained by sorting the height of all vegetation points relative to the ground in ascending order, i is the number of vegetation points in the DH data (1, 2, . . . , n), j is the percent. h mean , h min , and h max are the mean, minimum, and maximum canopy heights, respectively. CRR ranges from 0 to 1.

Accuracy Validation
We compared the accuracy of the forest canopy heights (CHM ATL08 , CHM ATL03 initial , CHM CPTC CHM CCR-1 , and CHM CCR-2 ) extracted from the ICESat-2 photon cloud data with that obtained from the CAF-LiCHy data. It has been shown that the vegetation height based on photon data is highly correlated with the 95th percentile height of conventional airborne LiDAR data [52]. Therefore, the 95th percentile height of the airborne CHM was used as the validation data to validate the accuracy of the data at the photon level and segment level. The determination coefficient (R 2 ), bias of mean (Bias), mean absolute error (MAE), and root mean square error (RMSE) were used to verify the accuracy of the canopy height obtained from ICESat-2 data. R 2 : Determination coefficient Bias: MAE: Mean absolute error MAE': Normalized mean absolute error RMSE: Root Mean Square Error where h ALS is the height obtained from the CAF-LiCHy data, h ATL is the height obtained by inversion from the ICESat-2 data, and n is the number of total photons. MAE' is a linear transformation of the MAE to obtain values in the range of [0, 1]. Figure 9 shows the canopy heights from ATL08, initial ATL03, and after cross-track photon correction and canopy height correction. The canopy height provided by ATL08 is significantly overestimated in this region compared with that extracted from the CAF-LiCHy LiDAR data. The MAE and RMSE exceed 20 m. Figure 9 and Table 3 show that although ATL08 provides the canopy height per 100 m segment, there is a substantial difference between its canopy height and that obtained from the airborne LiDAR, resulting in low accuracy. These results suggest that ATL08 may not be suitable for obtaining canopy heights in subtropical forest areas. The accuracy of CHMATL03 initial extracted by the conventional method from the ATL03 data is higher than that of CHMATL08, but the error is relatively large. After applying the cross-track photon correction before photon denoising and slope filtering for ground photon selection, the canopy height R 2 increased by 38.24%, MAE decreased by 35.71%, and the RMSE decreased by 35.98%, indicating a significant improvement in accuracy. The accuracy of canopy height estimation was further improved by applying the slope factor correction method. The MAE of CHMCCR-1 decreased by 15.58%, and the RMSE decreased by 12.73%. The MAE of CHMCCR-2 decreased by 23%, and the RMSE decreased by 19.23%. CHMCCR-2 showed better results than CHMCCR-1. The difference between CHMCCR-1 and CHMCCR-2 was due to the different Dd extracted from different ancillary images. As shown in Figure 9 and Figure 10 and Table 3, the canopy height correction further reduced the canopy height error. Compared to the CHMATL03 initial, the CHMCCR-2 improved the R 2 by 91.18% and reduced the MAE by 50.50% and the RMSE by 48.29%. The R 2 of CHMCCR-2 explained 65% of the variation in canopy height, indicating a significant improvement over the pre-correction estimates. The results of CHMCCR-1 and CHMCCR-2 show that the slope-based canopy correction method improves the goodness of fit and reduces the error.  The accuracy of CHM ATL03 initial extracted by the conventional method from the ATL03 data is higher than that of CHM ATL08 , but the error is relatively large. After applying the crosstrack photon correction before photon denoising and slope filtering for ground photon selection, the canopy height R 2 increased by 38.24%, MAE decreased by 35.71%, and the RMSE decreased by 35.98%, indicating a significant improvement in accuracy. The accuracy of canopy height estimation was further improved by applying the slope factor correction method. The MAE of CHM CCR-1 decreased by 15.58%, and the RMSE decreased by 12.73%. The MAE of CHM CCR-2 decreased by 23%, and the RMSE decreased by 19.23%. CHM CCR-2 showed better results than CHM CCR-1 . The difference between CHM CCR-1 and CHM CCR-2 was due to the different D d extracted from different ancillary images. As shown in Figures 9 and 10 and Table 3, the canopy height correction further reduced the canopy height error. Compared to the CHM ATL03 initial , the CHM CCR-2 improved the R 2 by 91.18% and reduced the MAE by 50.50% and the RMSE by 48.29%. The R 2 of CHM CCR-2 explained 65% of the variation in canopy height, indicating a significant improvement over the pre-correction estimates. The results of CHM CCR-1 and CHM CCR-2 show that the slope-based canopy correction method improves the goodness of fit and reduces the error.

Accuracy of Canopy Height Estimation for Different Segment Sizes and Relative Heights
This section describes the results of a sensitivity analysis of the forest canopy height parameters for different segment sizes and different percentile heights. Table 4 lists the canopy height accuracies before and after applying the proposed canopy height correction model for different segment sizes. The canopy height error decreases with an increase in the segment size, reaching the highest R 2 and the lowest MAE and RMSE for the 100 m segment. The larger the segment distance, the better the performance of the canopy height correction model is. As the segment size increases, the R 2 of CHMATL03 initial increases from 0. 35

Accuracy of Canopy Height Estimation for Different Segment Sizes and Relative Heights
This section describes the results of a sensitivity analysis of the forest canopy height parameters for different segment sizes and different percentile heights. Table 4 lists the canopy height accuracies before and after applying the proposed canopy height correction model for different segment sizes. The canopy height error decreases with an increase in the segment size, reaching the highest R 2 and the lowest MAE and RMSE for the 100 m segment. The larger the segment distance, the better the performance of the canopy height correction model is. As the segment size increases, the R 2 of CHM ATL03 initial increases from The sensitivity analysis of the forest parameters shows that the MAE and RMSE of CHM ATL03 initial show a rapid decrease with an increase in the segment size until 60 m, followed by a small decrease. The MAE and RMSE of CHM CTPC show a rapidly decreasing trend before the 50 m segment and a small decrease afterward. The MAE and RMSE of CHM CCR-2 show a rapidly decreasing trend before the 40 m segment and a small decrease afterward. The proposed cross-track photon correction and canopy height correction method improve the accuracy of canopy height estimation in a small area. In addition, CHM ATL03 initial , CHM CTPC , and CHM CCR-2 provide better performances than CHM ATL08 for the 100 m segment size. The percentile canopy height is crucial for describing the vertical structure of forests and modeling forest biomass. Therefore, an accuracy assessment of the percentile height was performed. Table 5 shows the sensitivity of eight percentiles of the forest canopy height at the photon level (RH70, RH75, RH80, RH85, RH90, RH95, RH98, and RH100). The R 2 range of CHM ATL08 is 0.17~0.19, that of CHM ATL03 initial is 0.33~0.36, that of CHM CTPC is 0.45~0.48, and that of CHM CCR-2 is 0.63~0.68. The R 2 does not change significantly with an increase in the percentile and is not shown in Table 5. The Bias, MAE, and RMSE of CHM ATL08 decrease with a decrease in the percentile canopy height. CHM ATL03 initial , CHM CTPC , and CHM CCR-2 exhibit a decrease followed by an increase in the Bias, MAE, and RMSE with the decreasing percentile canopy height. The difference is that CHM ATL03 initial reaches the minimum error at RH75, CHM CTPC at RH80, and CHM CCR-2 at RH98. The results show that the canopy-corrected heights of CHM CCR-2 have the lowest errors, indicating that the proposed canopy correction method reduces the error.

Effect of Slope on Canopy Correction
ICESat-2 was designed to overcome the disadvantages of ICESat-1 regarding the measurement error on steep slopes. However, the canopy heights estimated by ICESat-2 on different slopes reveal that the error increases with the slope. The Chinese forestry slope grading standard was used to classify the slope into six classes to quantify the influence of slope on canopy height estimation. As shown in Table 6, the R 2 of CHM ATL03 initial , CHM CTPC , and CHM CCR-2 decreases as the slope increases. The MAE of each slope class for CHM ATL03 initial is stable at 8.41~9.94 m, and the RMSE is stable at 11.50~13.10 m, showing no changes. The MAE and RMSE of CHM CTPC and CHM CCR-2 increase with the slope. The MAE of each slope class for CHM CTPC is stable at 2.94~4.24 m, and the RMSE is stable at 3.67~5.17 m. The MAE of each slope class for CHM CCR-2 is stable at 2.76~3.48 m, and the RMSE is stable at 3.37~4.41 m. The MAE and RMSE of CHM CCR-2 are on average 0.75 m and 0.84 m lower than that of CHM CTPC , respectively. The errors were lower after the cross-track photon correction and canopy height correction, indicating that these proposed methods reduced the canopy height errors due to slope. The effect of the cross-track photon correction was greater than that of the canopy height correction for reducing the error.

Comparison of ICESat-2 Canopy Height Inversion Accuracy
We proposed a canopy height correction method based on the cross-track photon correction of Li et al. [47] to minimize the error caused by the slope. This correction method may not be applicable to inclined trees, whose height measurements are more complicated [60]. Figure 9a shows that the ATL08 data significantly overestimate the canopy height in this region, which is consistent with the results of Sun et al. [61] but different from those of Narine et al. [21] and Duncanson et al. [53], who found that ATL08 data underestimated canopy heights. The likely reason is that ICESat-2 uses the Landsat Vegetation Continuous Fields (VCF) product [62] to define the forest. If the average VCF canopy cover per 10 km segment is >5%, the ATL08 algorithm assumes that this segment is forested [29]. Because land-use types have changed over the past decade, ATL08 may be using an algorithm for non-forested areas, resulting in anomalous height changes. In addition, due to the properties of the photon-counting LiDAR sensor and the ATL08 algorithm, it is possible that noisy photons were incorrectly identified as canopy or terrain photons in the 100 m segment size. This error would lead to an incorrect estimation of canopy height, especially for densely or sparsely vegetated areas [24]. Therefore, using canopy height data from ATL08 is not recommended in subtropical forest areas.
Using the 20 m segment size for photon denoising and photon classification significantly improved the canopy height estimates (Figure 9a). Streutker et al. [1] also showed that a small segment size is preferable for canopy photon selection. However, a smaller segment size is not necessarily better. A lower density of the photon cloud of the vegetation signal in small segments may classify the vegetation signal as noise photons. Compared to CHM CTPC and CHM CCR-2 , the canopy height of CHM ATL03 initial was significantly overestimated, probably due to noise near the ground or at the top of the canopy after photon denoising. This noise was removed by cross-track photon correction and SFMB filtering. CHM CCR-2 is based on the cross-track photon correction results and introduces a slope factor to correct the TOC photons. Compared with CHM CTPC , the MAE of CHM CCR-2 is reduced by 0.87 m and the RMSE is reduced by 0.90 m. The average correction value of TOC photons is 0.66 m. The reduced error of CHM CCR-2 is close to the average correction value of TOC photons. In areas with slope > 15 • , TOC photon correction reduces MAE by 0.76~0.97 m and RMSE by 0.76~1.09 m. R 2 improves by 0.18 after TOC photon correction, it reflects that the canopy correction method proposed in this paper can reduce the canopy height estimation error caused by slope. CHM CTPC without TOC photon correction had the lowest error at RH80. CHM CCR-2 with TOC photon correction had the lowest forest canopy height MAE and RMSE at RH98. It indicates that the TOC photon correction method eliminates the error caused by slope and restored the natural morphology of the forest canopy.
The ATL08 algorithm is suitable for global data processing, resulting in high errors in some regions. However, the canopy height estimates were improved after cross-track photon correction and canopy height correction based on ATL03 data, indicating that ICESat-2 could perform canopy height inversion in subtropical regions.

Influencing Factors of ICESat-2 CHM Inversion
This section describes the effects of cross-track photon correction, slope, segment size, percentile canopy height, D d , canopy density, canopy relief rate, growth, and harvesting on canopy height estimation using ICESat-2 data.
The effect of slope on the ICESat-2 canopy height estimation is reflected in two aspects. First, if the 3D photon data are converted into 2D photon profiles without cross-track photon correction [47], errors will occur in the DTM and DSM due to the slope influence; this problem is called the horizontal error. Second, the position of the trunk bottom may be overestimated or underestimated in the vertical direction (cf. Figure 6), which is called the vertical error. Li et al. [47] reported that the MAE and RMSE of the height value in the DTM were reduced after cross-track photon correction by 1.86~2.52 m and 2.12~2.50 m, respectively. The slope can affect the canopy height estimates obtained from the DTM, resulting estimates inaccurate. The relative canopy height is the height of the canopy photons minus the height of the fitted ground. Thus, on high slopes, the errors of the ground points are transferred to the canopy height estimates. As shown in Table 3, the MAE and RMSE of the canopy heights after cross-track photon correction were reduced by 2.15 m and 2.63 m, respectively, which is similar to the improvement in the DTM accuracy by cross-track photon correction. These results show that the cross-track photon correction method not only effectively improves the ground estimation accuracy but also the canopy height estimation accuracy. After the TOC photon correction, the MAE was further reduced by 0.89 m, and the RMSE was further reduced by 0.90 m, suggesting that the horizontal error is larger than the vertical error and that accurate terrain estimation is crucial for deriving canopy heights.
In mountainous regions, terrain undulations cause substantial uncertainty in canopy height estimation [63]. We quantified the error in the vertical direction caused by the slope. Figure 11a shows the influence of the slope on the absolute correction value. On south-facing slopes, the difference was negative when the TOC photons were north of the canopy centroid and positive when they were south of the canopy centroid. The opposite is true for north-facing slopes. The difference between the CHMs before and after correction ranged from −2.55 m to 2.90 m. A tree canopy is typically symmetrical with the trunk as the axis, and the positive and negative correction values are determined by the relative positions. Here, the correction values are the absolute values to depict the effect of slope on canopy correction. As the slope becomes steeper, the absolute correction value increases. On slopes ≤ 10 • , the absolute correction values are less than 0.50 m, indicating that the canopy height correction method is not required in low-slope areas. As shown in Table 6, the canopy heights estimated after cross-track photon correction had relatively high accuracy in flat areas. The absolute correction values are different for the same slope values due to the difference in D d . Moreover, the fitted curve in Figure 11a shows that the vertical error increases exponentially with the slope. Therefore, the effect of the topographic slope must be considered in very steep areas. Future studies should evaluate the empirical equation of the slope correction method for ICESat-2 canopy height estimate under different environments. south-facing slopes, the difference was negative when the TOC photons were north of the canopy centroid and positive when they were south of the canopy centroid. The opposite is true for north-facing slopes. The difference between the CHMs before and after correction ranged from −2.55 m to 2.90 m. A tree canopy is typically symmetrical with the trunk as the axis, and the positive and negative correction values are determined by the relative positions. Here, the correction values are the absolute values to depict the effect of slope on canopy correction. As the slope becomes steeper, the absolute correction value increases. On slopes ≤ 10°, the absolute correction values are less than 0.50 m, indicating that the canopy height correction method is not required in low-slope areas. As shown in Table 6, the canopy heights estimated after cross-track photon correction had relatively high accuracy in flat areas. The absolute correction values are different for the same slope values due to the difference in Dd. Moreover, the fitted curve in Figure 11a shows that the vertical error increases exponentially with the slope. Therefore, the effect of the topographic slope must be considered in very steep areas. Future studies should evaluate the empirical equation of the slope correction method for ICESat-2 canopy height estimate under different environments. It is necessary to meet certain accuracy requirements for smaller segment sizes considering the future demand for forest parameters. We compared the sensitivity of the forest parameters for 10 segment sizes (10 m to 100 m) (cf. Table 4). The accuracy increases with an increase in the segment length, in agreement with the findings of Streutker et al. It is necessary to meet certain accuracy requirements for smaller segment sizes considering the future demand for forest parameters. We compared the sensitivity of the forest parameters for 10 segment sizes (10 m to 100 m) (cf. Table 4). The accuracy increases with an increase in the segment length, in agreement with the findings of Streutker et al. [1]. It is expected that the MAE and RMSE will continue to decrease as the segment size increases but at a slower rate. The R 2 was 0.64, and the RMSE was 3.64 m for the 30 m segment size after cross-track photon correction. Our R 2 and RMSE values were slightly smaller than those of Zhu et al. (0.70 and 4.31 m, respectively [64]. After the cross-track photon correction and canopy height correction, the R 2 was 0.70, and the RMSE was 3.21 m, better than the results of Zhu et al. [64]. The experimental results show that the combination of cross-track photon correction and canopy height correction substantially reduces the canopy height estimation error due to slope for ICESat-2 data. A comparison of the different percentile canopy heights showed that the corrected values were more accurate than the uncorrected values, and the MAE and RMSE were the lowest at RH98. The forest canopy height results without cross-track photon correction and TOC photon correction were optimal at RH75, consistent with the findings of Neuenschwander and Pitts [29]. The forest canopy height results with cross-track photon correction were optimal at RH80. This finding indicates that canopy height correction can reduce the uncertainty in selecting the optimal percentile canopy height for different terrain conditions. Figure 11b shows the proportions of the intervals for D d = 0.50 m and the maximum, average, and minimum corrected heights. The TOC photon selection algorithm and the crown segmentation algorithm affect D d . The highest point or percentile methods are typically used for TOC photon selection [19,20,26,28,33], and the obtained TOC photon results are relatively fixed. The canopy segmentation results depend on the tree species, canopy width, canopy morphology, planting density, solar position, and shading. The segmentation procedure is dependent on the user-specified parameters (scale, shape ratio, and compactness ratio). When trees are grouped closely together in dense stands, the ends of their branches may touch the crowns of nearby trees, resulting in larger crowns [65]. Thus, the crown centroids may have errors, which are transferred to D d , leading to uncertainty in canopy height correction. This problem causes the largest uncertainty in the proposed canopy correction method. It is expected that crown segmentation algorithms will be improved using deep learning and other methods, minimizing this problem. Figure 11b shows that the proportion of D d < 2.0 m accounts for 76.96%, which is attributed to the canopy width of 4.0~8.0 m in this area. The maximum correction value does not occur in the interval with the maximum D d . Equation (7) indicates that the error is caused by the slope. This result suggests that the effect of slope on canopy height is not negligible, and the larger the slope, the larger the D d is, and the larger the canopy height error is. Figure 11c shows the normalized MAE of CHM CCR-2 at different canopy densities. The dashed line is the boundary line of normalized MAE at different canopy densities. Zhu et al. [34] also found that forest canopy density substantially affected the inversion of canopy height from ICESat-2 data. We classified the canopy density into seven classes to investigate the canopy density effect on the canopy height estimates of ICESat-2 data. The mean normalized MAE increased with increasing canopy density, in agreement with the results of Narine [21], Neuenschwander and Magruder [24], and Martin [66]. This study differs from others in that the MAE was not always high in areas with high canopy density, but the uncertainty of the normalized MAE rose with the increasing canopy density. The boundary line shows the maximum error of the seven canopy density classes, indicating an increasing trend of the normalized MAE. Although there is some uncertainty in this line, it highlights the differences in the error for each canopy density class, an approach rarely used in data analysis.
Existing studies on extracting forest vegetation parameters from photon-counting LiDAR data have mostly focused on canopy height inversion. Canopy height is a crucial parameter. Our study with ICESat-2 data only considered the canopy segments but not the internal canopy characteristics, and a CHM cannot describe the canopy geometry. Therefore, we used the CRR to analyze the geometric characteristics of the tree crowns to obtain a macroscopic description of the canopy height characteristics. Because CHM CCR-2 provided better performance than CHM CCR-1 , this section compares the CRR of CHM CCR-2 and CHM ALS . The CRR of the spaceborne LiDAR data was significantly correlated with the airborne LiDAR data, with a Pearson correlation coefficient of 0.88 and an R 2 of 0.77 (cf. Figure 11d). This result shows that the CHM derived from the spaceborne LiDAR data and the airborne LiDAR data are similar in terms of canopy characteristics. The CRR is not related to the accuracy of the ICESat-2 inversion of canopy height but indicates that the ICESat-2 CHM is relatively accurate. The CRR can be used as an index to demonstrate the accuracy of ICESat-2 data for characterizing the canopy in addition to other accuracy metrics (MAE, RMSE, etc.).
There was a difference in the acquisition time between the spaceborne LiDAR data and the airborne LiDAR data. The airborne LiDAR data were collected in April 2018, and the ICESat-2 data were collected in October 2019; thus, an exact match in time was impossible. During the one-and-a-half-year period, the trees grew, and their height changed, leading to errors. In addition, this area is an operational forest, where logging operations occur, resulting in height estimation errors. Tree planting in non-forested areas can also cause errors in canopy height estimation. Vegetation planting, growth, and harvesting affect the accuracy of forest canopy height estimation.

Future Directions and Implications
Forest height indicators extracted from ICESat-2 data have been used as critical variables for estimating biomass [21,53,67]. Our results suggest that using ICESat-2 canopy height indicators is more challenging in stands with complex terrain and high variability in canopy density. In contrast, cross-track photon correction and canopy height correction can improve the accuracy of canopy height estimates.
The findings of this study show that ICESat-2 data can be used successfully to assess the structure and characteristics of forest stands with different canopy heights, enabling large-scale monitoring and management. The results obtained from the Google Earth satellite images show the potential of their use for validation in different regions. Martin et al. [68] detected more than 1.8 billion trees in the Sahara, Sahel, and sub-humid zone of West Africa using high-resolution satellite images and depth learning algorithms. Although this method was applied to savanna areas, forest canopy segmentation on a global scale will become a reality with the development of ultra-high spatial resolution satellites and single-tree identification algorithms. Canopy correction methods will be further improved. Any single detection system has some limitations. Thus, combining different ultra-high spatial resolution satellite images (Worldview-2/3, GeoEye-1, etc.) with ICESat-2 data can be advantageous. This study demonstrates the capacity to utilize these data for estimating forest attribute values and provides a thorough review of extracting forest structure parameters from ICESat-2 data. Nie et al. [69] developed a theoretical model to quantify the effects of terrain slope, canopy radius, canopy shape, and offset distance to the slope on canopy height displacement in ellipsoidal canopies. The effect of slope on canopy height estimation may differ for different types of canopies and terrain. This paper proposes a theoretical approach to eliminating slope errors. In the future, we can build on this theory to verify the effect of different canopy shapes on canopy height estimation experimentally. ICESat-2 will be integrated with other datasets to extract detailed estimates of forest inventory parameters for monitoring purposes on a continuous time scale. This strategy is applicable due to the growing availability of ICESat-2 data and the improved capabilities for forest monitoring.

Conclusions
Canopy height inversion from ICESat-2 laser altimetry data is substantially influenced by slope, canopy density, and signal-to-noise ratio. Studies have shown that slope is the most influential factor. Our results indicate that the vertical error increases exponentially with the increasing slope; thus, the effect of slope on canopy height estimation must be considered. This study proposed the first ICESat-2 CCR based on high spatial resolution DOM. The effect of slope on the canopy height estimates obtained from ICESat-2 data was minimized by the cross-track photon correction method. This method reduced the MAE by 35.71% and the RMSE by 35.98%. The canopy correction method based on high spatial resolution DOM based on the cross-track photon correction resulted in a further 23% reduction in the MAE and a 19.23% reduction in the RMSE. The proposed method significantly improved the accuracy of ICESat-2 forest canopy height estimation compared with the traditional method. Using uncorrected canopy height data from ATL08 is not recommended in subtropical forest areas.
This study investigated the sources of slope error and their effect on canopy height estimation and proposed a slope error correction method. The method is suitable for spaceborne photon-counting LiDAR systems and improves the accuracy of canopy height estimation. However, the correction method uses multi-source data, increasing the computational complexity. Thus, the method is more time-consuming than using ATL08 data for large photon point clouds. Since we focused on a theoretical approach to minimize slope errors, we did not quantify the effect of canopy morphology on canopy height estimation. Future studies will focus on the validation of our theoretical approach in different regions and ecosystems to demonstrate the applicability of the methodology.