Mapping the Forest Canopy Height in Northern China by Synergizing ICESat-2 with Sentinel-2 Using a Stacking Algorithm

: The forest canopy height (FCH) plays a critical role in forest quality evaluation and resource management. The accurate and rapid estimation and mapping of the regional forest canopy height is crucial for understanding vegetation growth processes and the internal structure of the ecosystem. A stacking algorithm consisting of multiple linear regression (MLR), support vector machine (SVM), k -nearest neighbor ( k NN), and random forest (RF) was used in this paper and demonstrated optimal performance in predicting the forest canopy height by synergizing Sentinel-2 images acquired from the cloud-based computation platform Google Earth Engine (GEE) with data from ICESat-2 (Ice, Cloud, and Land Elevation Satellite-2). This research was conducted to achieve continuous mapping of the canopy height of plantations in Saihanba Mechanical Forest Plantation, which is located in Chengde City, northern Hebei province, China. The results show that stacking achieved the best prediction accuracy for the forest canopy height, with an R 2 of 0.77 and a root mean square error (RMSE) of 1.96 m. Compared with MLR, SVM, k NN, and RF, the RMSE obtained by stacking was reduced by 25.2%, 24.9%, 22.8%, and 18.7%, respectively. Since Sentinel-2 images and ICESat-2 data are publicly available, this opens the door for the accurate mapping of the continuous distribution of the forest canopy height globally in the future.


Introduction
The forest canopy height plays a significant role in vegetation health assessment [1,2], and is also critical for deriving several key biophysical parameters, such as the aboveground biomass (AGB), carbon storage, biodiversity, vegetation photosynthesis, carbon cycle, and global climate change [3][4][5]. Obtaining accurate, large-scale, and multitemporal forest canopy height information can provide policymakers with scientific information on deforestation, forest degradation, and forest quality [6].
The application of remote sensing technology in forest surveying has greatly improved the efficiency of forest management [7]. Optical remote sensing, synthetic aperture radar (SAR), and light detection and ranging (LiDAR) data have been widely used in land change and vegetation cover mapping [8][9][10][11][12][13]. Among these, LiDAR has proved to be the most accurate in estimating forest vertical structure parameters, with its unique advantage combine multiple base models that are functionally independent classifiers or regressors to obtain a better and more comprehensive prediction model [7,38]. Even if one base learner gives a wrong prediction, other base models can correct and improve it [42]. As the representative of ensemble learning, the stacking algorithm is able to train a model to combine other base models. Stacking can integrate the advantages of each individual base model to establish a prediction model with a better estimation effect and applicability [43]. Li et al. [44] mapped the growing stem volume (GSV) of plantations in northeast China successfully using the stacking approach and achieved satisfactory estimation results. In addition, the algorithm has also been applied to predict soil organic matter [45].
This study aimed to develop a stacking approach for estimating and mapping the forest canopy height in northern China with satisfactory accuracy. Specifically, this study combined ICESat-2 data with synthetic Sentinel-2 images obtained from GEE for developing of the stacking, and the derived forest heights were validated with ground-truth data collected in Saihanba, Heibei, China.

Study Area
The Saihanba Mechanical Forest Plantation is located in Chengde City, northern Hebei province, China (116 • 51 ~117 • 39 E, 42 • 02 ~42 • 36 N) ( Figure 1) and covers an area of 933.33 km 2 , with an elevation ranging from 1010 to 1940 m. The study area is characterized by a temperate continental monsoon climate, with a minimum temperature of −43.2 • C and an average annual snow cover of 7 months. The annual mean precipitation is about 490 mm, with an average annual frost-free period of 68 days. The Saihanba forest farm was established in 1962, and the total volume of trees has reached 10.12 million m 3 currently. The dominant tree species are larch (Larix ologensis), Scots pine (Pinus sylvestris), birch (Betula platyphylla Suk), and spruce (Picea asperata Mast). used in estimating vegetation parameters [34][35][36]. Compared with individual models, ensemble learning is able to combine multiple base models that are functionally independent classifiers or regressors to obtain a better and more comprehensive prediction model [7,38]. Even if one base learner gives a wrong prediction, other base models can correct and improve it [42]. As the representative of ensemble learning, the stacking algorithm is able to train a model to combine other base models. Stacking can integrate the advantages of each individual base model to establish a prediction model with a better estimation effect and applicability [43]. Li et al. [44] mapped the growing stem volume (GSV) of plantations in northeast China successfully using the stacking approach and achieved satisfactory estimation results. In addition, the algorithm has also been applied to predict soil organic matter [45].
This study aimed to develop a stacking approach for estimating and mapping the forest canopy height in northern China with satisfactory accuracy. Specifically, this study combined ICESat-2 data with synthetic Sentinel-2 images obtained from GEE for developing of the stacking, and the derived forest heights were validated with ground-truth data collected in Saihanba, Heibei, China.

Study Area
The Saihanba Mechanical Forest Plantation is located in Chengde City, northern Hebei province, China (116°51′~117°39′ E, 42°02′~42°36′ N) ( Figure 1) and covers an area of 933.33 km 2 , with an elevation ranging from 1010 to 1940 m. The study area is characterized by a temperate continental monsoon climate, with a minimum temperature of −43.2 °C and an average annual snow cover of 7 months. The annual mean precipitation is about 490 mm, with an average annual frost-free period of 68 days. The Saihanba forest farm was established in 1962, and the total volume of trees has reached 10.12 million m 3 currently. The dominant tree species are larch (Larix ologensis), Scots pine (Pinus sylvestris), birch (Betula platyphylla Suk), and spruce (Picea asperata Mast).

Measurement of the Forest Canopy Height in Saihanba
The forest management inventory (FMI), which is carried out every ten years in China, plays a key role in forest resources monitoring and evaluation. A forest resource

Measurement of the Forest Canopy Height in Saihanba
The forest management inventory (FMI), which is carried out every ten years in China, plays a key role in forest resources monitoring and evaluation. A forest resource database Remote Sens. 2021, 13, 1535 4 of 17 composed of irregular blocks was established and is updated once a year in the FMI to guide the management of forests [26,46]. The database was updated in 2017, including land type, block area, forest canopy height, and tree species, and this information was analyzed to obtain the forest canopy information of Saihanba.
In our study, 5404 sample plots composed of seven tree species were selected from Saihanba forest experiment sites ( Figure 2). As the most widely distributed tree species, larch accounted for 56.2% of all sample plots. Its canopy height ranged from 3.7 m to 20.0 m, and the mean, standard deviation, and coefficient of variation were 15.39 m, 3.69 m, and 23.9%, respectively. The canopy height of spruce had the largest coefficient of variation of all the tree species, while the mean value and standard deviation of oak was the lowest ( Table 1). database composed of irregular blocks was established and is updated once a year in the FMI to guide the management of forests [26,46]. The database was updated in 2017, including land type, block area, forest canopy height, and tree species, and this information was analyzed to obtain the forest canopy information of Saihanba. In our study, 5404 sample plots composed of seven tree species were selected from Saihanba forest experiment sites ( Figure 2). As the most widely distributed tree species, larch accounted for 56.2% of all sample plots. Its canopy height ranged from 3.7 m to 20.0 m, and the mean, standard deviation, and coefficient of variation were 15.39 m, 3.69 m, and 23.9%, respectively. The canopy height of spruce had the largest coefficient of variation of all the tree species, while the mean value and standard deviation of oak was the lowest (Table 1).

ICESat-2 ATLAS Data and Processing
The ICESat-2 carries an advanced terrain laser altimeter system (ATLAS), which emits three pairs of beams. The distance between each pair of beams is about 3.3 km, and the distance between two trajectories in the beam is 90 m. Three pairs of laser beams represent six ground trajectories, and the corresponding data numbers from left to right are, respectively, Ground Track 1 Left (GT1L) to Ground Track 3 Right (GT3R) [14,19,47,48].
The ICESat-2 products include three levels, from ATL01 to ATL21, and the ATL08 data (level 3) covering the whole study area from 5 November 2018 to 8 December 2019 were downloaded from the National Ice and Snow Data Center (NSIDC) (https://nsidc.org/ data/ATL08/ accessed on 8 March 2021) (Figure 3a). The ATL08 product directly provides data on the terrain, height metrics, and apparent surface reflectance (ASR) (Table 2) in fixed 100 m data segments along the track, with a diameter of 17 m ( Figure 3) [23]. The height metrics include minimum height, mean height, median height, maximum height, and 10 height profile quantiles, which can be used as variables for forest height modeling and estimation, together with ASR.     Terrain, height metrics, and apparent surface reflectance attributes from the ATL08 product for our study area were extracted. The canopy height was calculated from 98% of all individual relative canopy heights in the region. As a PCL system, ATL08 is easily affected by the noise from the atmosphere and the sun background. Therefore, ATL08 was denoised before processing and obtaining the vegetation products [49]. Segments of less than 50 signal photons were filtered out, as they did not characterize the ground surface accurately. To make full use of available ATL08, the differential, regressive, and Gaussian adaptive nearest neighbor (DRAGANN) methods were developed to identify and remove noise photons from the histogram of photon point clouds [49,50]. In addition, the segments with an average canopy height greater than 50 m and less than 2 m were excluded as invalid values. Table 3 shows the canopy height distribution for all 5404 segments covering the study area. The canopy height of the segments ranged from 15 to 20 m, accounting for 48.5% of all segments. By contrast, the canopy height for all the segments varied broadly from 4.26 to 26.41 m, with a mean, standard deviation, and coefficient of variation of 15.11 m, 3.62 m, and 23.9%, respectively.

Sentinel-2 Image Preprocessing and Spectral Variable Calculation
While ICESat-2 cannot provide wall-to-wall coverage, combining it with continuous remote sensing data, such as those from Sentinel-2, proved to be an effective method with which to map the forest canopy heights [17]. Sentinel-2 was launched in June 2015, and it is now widely used in forest parameter estimation and land cover mapping, as the red-edge bands are sensitive to the chlorophyll of vegetation leaves [23]. Level-2A Sentinel-2 images in the growing season (June-September) from 2018 to 2019 were downloaded with the cloud-based computation platform Google Earth Engine (GEE, https://code.earthengine. google.com/ accessed on 8 March 2021). The Sentinel-2 images obtained were a Level-2A product; they were atmospherically corrected by the Sen2Cor module, and the pixels were transformed to account for surface reflectance. The quality assessment band (QA60) was used for cloud masking [26]. To obtain more stable and higher resolution pixel values, the visible bands (red, green, blue, and near-infrared), with a resolution of 10 m, and three red-edge bands (RedEdge 1-3), with a resolution of 20 m, were selected and applied to obtain a median image composition. A total of eight vegetation indices, including five common vegetation indices and three red-edge vegetation indices, were extracted as spectral features to establish a link between ICESat-2 and the mapping (Table 4).

Vegetation Index Description Reference
Normalized difference vegetation index (NDVI) 2.5. Model Development and Assessment 2.5.1. Parametric and Nonparametric Models Parametric and nonparametric models are commonly used for forest parameter estimation [30][31][32]35]. Parametric models, such as multiple linear regression (MLR), are easy to implement and can describe the quantitative relationship between multiple variables and the forest canopy height. However, under complex forest conditions, the relationship between remote sensing variables and the forest canopy height may be complex [32]. Nonparametric models have the potential to estimate forest parameters with satisfactory accuracy and have been proved to be effective [35]. This study used support vector machine (SVM), k-nearest neighbor (kNN), and random forest (RF) to estimate forest height compared with MLR. In addition, a stacking algorithm composed of these four models was established to explore its potential in estimating the forest canopy height.
MLR assumes that there is a highly linear relationship between the dependent variable and independent variables; it can measure the effect of multiple variables on the prediction, and it is a commonly used parametric model [30]. It is more effective and practical to predict or estimate the dependent variable with a combination of multiple independent variables than by using only one independent variable [31]. The general expression of a multiple linear regression equation is shown in Equation (1).
where Y is the dependent variable, X 1 , X 2 , X 3 . . . X n is the independent variable, β 1 , β 2 , β 3 β n is the regression coefficient, and C is a constant. SVM can achieve good prediction performance, even if there are few samples available. The nonlinear kernel function can transform the input data into high-dimensional feature space and reduce the error and complexity of the model [36,40,41]. The radial basis function, Gaussian function, and exponential kernel function were established and compared, and the kernel function with the minimum estimation error was determined. The "tune.svm" function in R (package "e1071" [55]) can be used to set the initial values of the penalty coefficient (cost) and kernel function parameter (gamma) (usually 10 −3 -10 3 ), and then the value with the minimum error rate is selected as the optimal parameter [34].
The kNN algorithm searches the k samples nearest to the sample plot to be predicted by spectral distance, and then uses the specific prediction rules and the attributes of the selected samples to predict the sample plot to be predicted [37]. Samples with a normal distribution are not required in kNN, making it flexible, and it is widely used in estimating forest parameters and biodiversity mapping. The distance measurement, k value, and prediction rules can strongly affect the estimation results [36]. In Saihanba, the Euclidean distance and reciprocal weighting of the distance method were used to predict the forest canopy height (Equation (2)). The k value was set from 1 to 50, and the number of samples with the minimum estimation error was determined. where FCH p denotes the predicted FCH at the pixel P, y j denotes the FCH observation corresponding to the jth sample plot, d pj is the spectral distance from the pixel P to the jth sample plot, and k denotes the optimal number of sample plots. RF is a widely used nonparametric algorithm in land classification and forest parameter estimation because it is robust, relatively easy to understand, and insensitive to noise data [38]. In addition, it is an ensemble learning method; it constructs a large number of regression trees for the prediction. The number of regression trees (ntree, with a default value of 500) and the number of nodes (mtry, within the number of input variables) are the basic parameters of random forests. The package "randomForest" [56] in R was used for the RF process.

The Stacking Algorithm
The stacking algorithm combines the advantages of the base models, which can greatly improve the prediction effect [44]. First, stacking used the base models to combine training samples to build its own prediction models; then, all the prediction results were set as independent variables to conduct new joint modeling. Stacking generalizes the output of all base models, which can significantly improve the prediction accuracy and reduce the overall error.
In Saihanba, a stacking algorithm consisting of MLR, SVM, kNN, and RF was developed to estimate and map the forest canopy height (Figure 4). The four base models were also used for a direct comparison. diction rules can strongly affect the estimation results [36]. In Saihanba, the Euclidean distance and reciprocal weighting of the distance method were used to predict the forest canopy height (Equation (2)). The k value was set from 1 to 50, and the number of samples with the minimum estimation error was determined.
where denotes the predicted FCH at the pixel P, y denotes the FCH observation corresponding to the jth sample plot, is the spectral distance from the pixel P to the jth sample plot, and k denotes the optimal number of sample plots.
RF is a widely used nonparametric algorithm in land classification and forest parameter estimation because it is robust, relatively easy to understand, and insensitive to noise data [38]. In addition, it is an ensemble learning method; it constructs a large number of regression trees for the prediction. The number of regression trees (ntree, with a default value of 500) and the number of nodes (mtry, within the number of input variables) are the basic parameters of random forests. The package "randomForest" [56] in R was used for the RF process.

The Stacking Algorithm
The stacking algorithm combines the advantages of the base models, which can greatly improve the prediction effect [44]. First, stacking used the base models to combine training samples to build its own prediction models; then, all the prediction results were set as independent variables to conduct new joint modeling. Stacking generalizes the output of all base models, which can significantly improve the prediction accuracy and reduce the overall error.
In Saihanba, a stacking algorithm consisting of MLR, SVM, kNN, and RF was developed to estimate and map the forest canopy height (Figure 4). The four base models were also used for a direct comparison.

Model Assessment
The ICESat-2 segments were randomly divided into two groups: 70% (training set, n = 3782) for model training and 30% (testing set, n = 1622) for independent validation. The coefficient of determination (R 2 ) and the RMSE [57] were calculated to evaluate the model performance with the testing set.
Remote Sens. 2021, 13, 1535 where y i is the observed value,ŷ i is the predicted value based on the models, y is the mean of all the observed values, and n is the sample size. All of the models and calculations were carried out using R 4.0.3 software.

Forest Canopy Height Estimation Models with ICESat-2
All the parameters were significantly and positively correlated with the forest canopy height (p < 0.001) ( Table 5). Among them, the 98th percentile height had the highest correlation coefficient of 0.71. The Pearson correlation coefficient was used to detect linear relationships between the variables and forest height. However, the relationship between the variables and height is not always simply linear in dense forest situations. To add more variable evaluation methods, importance evaluation based on random forest was added to explore the nonlinear relationship between ASR, height metrics, and the measured forest canopy height. The random forest algorithm was used to evaluate the importance of all variables, and the importance ranking was obtained, as shown in Figure 5. The 98th percentile height had the highest importance, while the minimum had the lowest.
Remote Sens. 2021, 13, x FOR PEER REVIEW 10 of 17 Figure 5. Importance ranking (with standard deviation) of the ASR and height metrics obtained using random forest. Figure 5. Importance ranking (with standard deviation) of the ASR and height metrics obtained using random forest.
The performances of MLR, SVM, kNN, and RF were similar. Stacking achieved the highest fitting and the lowest error, and the RMSE from the stacking model was reduced by 25.2%, 24.9%, 22.8%, and 18.7% compared with each individual model. Figure 6 shows the forest canopy height values predicted by the five models versus the observed values. There were different degrees of overestimation and underestimation when predicting the forest canopy height in Saihanba. However, compared with the other models, stacking greatly reduced overestimation and underestimation, achieving the best fit. The residuals obtained by stacking tended to be randomly distributed.

ICESat-2 Derived Estimation for Sentinel-2
The importance of variables obtained from the RF model showed that ARVI contributed the most to the model, followed by RECI and NDVI (Figure 7a). The least important variable in the model was the enhanced vegetation index (EVI). The variables were selected with the random forest algorithm based on importance ranking, and the final combination of variables was used to map the forest canopy height continuously, combined with the stacking algorithm. The spatial pattern of the forest canopy height predicted by

ICESat-2 Derived Estimation for Sentinel-2
The importance of variables obtained from the RF model showed that ARVI contributed the most to the model, followed by RECI and NDVI (Figure 7a). The least important variable in the model was the enhanced vegetation index (EVI). The variables were selected with the random forest algorithm based on importance ranking, and the final combination of variables was used to map the forest canopy height continuously, combined with the stacking algorithm. The spatial pattern of the forest canopy height predicted by stacking was highly consistent with those from ICESat-2 (Figure 7b). By extracting continuous optical features from Sentinel-2 and combining them with stacking, the ICESat-2 derived forest canopy height could be predicted with satisfactory accuracy. The R 2 and RMSE of the estimation model were 0.92 and 1.37 m, respectively (Figure 7b). stacking was highly consistent with those from ICESat-2 (Figure 7b). By extracting continuous optical features from Sentinel-2 and combining them with stacking, the ICESat-2 derived forest canopy height could be predicted with satisfactory accuracy. The R 2 and RMSE of the estimation model were 0.92 and 1.37 m, respectively (Figure 7b).

Discontinuous and Continuous Mapping
ICESat-2 can provide fixed spacing footprints for estimating the forest canopy height with satisfactory accuracy. However, discontinuous footprints are limited for predicting a comprehensive spatial distribution (Figure 8a). The Sentinel-2 images, which were obtained in the growing season of the vegetation, were used as the medium to extract continuous optical variables and were combined with stacking to map the forest canopy height distribution of Saihanba (Figure 8b).
Most of the larger forest canopy height prediction values were distributed in the northeast and south, while the values in the central region were relatively small. The minimum forest canopy height values were distributed in the western part of the region (Figure 8). The spatial pattern of the predicted forest canopy height followed the actual distribution in Saihanba.

Discontinuous and Continuous Mapping
ICESat-2 can provide fixed spacing footprints for estimating the forest canopy height with satisfactory accuracy. However, discontinuous footprints are limited for predicting a comprehensive spatial distribution (Figure 8a). The Sentinel-2 images, which were obtained in the growing season of the vegetation, were used as the medium to extract continuous optical variables and were combined with stacking to map the forest canopy height distribution of Saihanba (Figure 8b). Most of the larger forest canopy height prediction values were distributed in the northeast and south, while the values in the central region were relatively small. The minimum forest canopy height values were distributed in the western part of the region (Figure 8). The spatial pattern of the predicted forest canopy height followed the actual distribution in Saihanba.

ICESat-2 Metrics for Canopy Height Estimation
Height metrics extracted from ICESat-2 ATL08 are crucial for improving the efficiency of the stacking model and achieving high-precision canopy height estimation [17,33]. However, ICESat-2 can provide only a limited number of height metrics, and thus all variables need to be employed in the stacking model. Our results showed that the height metrics and ASR extracted from ICESat-2 were closely related to the forest canopy height. The lowest RMSE was achieved by using these variables together with the stacking model. Although the correlation between the ASR, minimum height, and forest canopy height were also significant (p < 0.001), the correlation coefficients were less than 0.1. To verify the effectiveness of these two variables for estimating the forest canopy height, the rest of the variables (excluding the ASR and minimum height) were used to established models to estimate the forest canopy height. The results showed that stacking still achieved the best fitting, and the smallest RMSE. Compared with MLR, SVM, kNN, and RF, the RMSE of stacking was reduced by 20.1%, 21.0%, 20.4%, and 18.3%, respectively. However, after excluding the ASR and minimum height, the RMSE of all models increased sharply. By contrast, the RMSE of stacking increased by 8.4%, indicating that the two variables significantly improved the estimation of the forest canopy height.
In addition, Narine et al. [14] used stepwise regression combined with the collinear elimination method to select variables from all height metrics for AGB estimation and obtained satisfactory accuracy. The method was also used in Saihanba, and the RMSE of the linear variable selection method was 10.9% higher than the results in Table 5, although the best model was still stacking. The height metrics provided by ATL08 were mainly percentile values of canopy height, which may lead to collinearity among variables. In complex forest conditions, the relationship between variables and forest parameters is not simply linear, so collinearity may not be crucial for nonparametric models with a higher estimation accuracy [31]. However, due to the limited number of variables and high correlation, modeling with all the variables extracted from ATL08 did not lead to redundancy and reduced the estimation efficiency.

Uncertainty, Limitations, and Suggestions
The accuracy of modeling of the forest canopy height is influenced by several factors, such as the topographic relief, tree species, and canopy cover [7,24,26]. Firstly, with the increase of the forest canopy cover, a limited number of photons can reach the ground and be reflected back to the sensor, thus limiting the application of ICESat-2 in forests with high canopy cover [18]. To reduce the uncertainty associated with high canopy cover, footprints with few photons were excluded from the forest canopy height modeling. Second, the forest canopy height model may vary from one species to another.
Out results showed that the poplar (n = 17), Chinese pine (n = 3), and spruce (n = 23) had a small number of samples in the testing set, and the RMSEs were 1.44, 0.45, and 2.41 m, respectively. Larch and birch accounted for the largest number of samples, with RMSEs of 1.96 and 1.95 m, respectively, which was reasonable given the overall estimation accuracy ( Figure 9). However, oak had the largest estimation error. Similar to spruce, the accuracy of oak was limited due to the large variation coefficients (33.6% and 26.4%) of the samples.
Out results showed that the poplar (n = 17), Chinese pine (n = 3), and spruce (n = 23) had a small number of samples in the testing set, and the RMSEs were 1.44, 0.45, and 2.41 m, respectively. Larch and birch accounted for the largest number of samples, with RMSEs of 1.96 and 1.95 m, respectively, which was reasonable given the overall estimation accuracy ( Figure 9). However, oak had the largest estimation error. Similar to spruce, the accuracy of oak was limited due to the large variation coefficients (33.6% and 26.4%) of the samples. It is also worth noting that ground-truth data were not collected in the same season and year as the ICESat-2 data. In our study, the ground-truth forest canopy height was measured in 2017, while the ICESat-2 data covering Saihanba was acquired between November 2018 and December 2019. This mismatch between ground-truth data collection It is also worth noting that ground-truth data were not collected in the same season and year as the ICESat-2 data. In our study, the ground-truth forest canopy height was measured in 2017, while the ICESat-2 data covering Saihanba was acquired between November 2018 and December 2019. This mismatch between ground-truth data collection and satellite data acquisition may have affected the forest canopy modeling.
The ICESat-2 data does not provide wall-to-wall spatial coverage [14]. Comprehensive mapping of the forest canopy height requires optical data for extrapolation [18]. As a powerful cloud platform, GEE can provide a variety of time series data, such as those from Landsat, Sentinel, and MODIS, to synthesize more stable median images after cloud masking [27]. Combining ICESat-2 with continuous coverage optical images has become an effective method to obtain a high-precision forest canopy height [18,19]. In addition, the Global Ecosystem Dynamics Investigation (GEDI) can provide vertical canopy information between 52 • N and S latitudes, which adds another source of LiDAR data [19,58]. Integrating ICESat-2, GEDI, and LiDAR with other forms of optical data, as well as environmental variables, and using machine learning algorithms to estimate the forest canopy height, may have significant potential in the future [59][60][61][62][63][64][65].

Conclusions
The accurate estimation of the forest canopy height is of great significance for assessing forest health and terrestrial ecosystem dynamics. In this study, a stacking algorithm consisting of MLR, SVM, kNN, and RF was developed to map the spatial pattern of the forest canopy height in the Saihanba forest experiment site in northern China by coupling the ICESat-2 and synthesized Sentinel-2 time series imagery. The results demonstrated that: (1) the canopy height parameters extracted from ICESat-2 were highly correlated with the field measurements. The correlation between the height metrics and the measured forest canopy height was higher than that of ASR, the 98th percentile height achieved the best correlation (R = 0.71), and (2) compared to individual algorithms (MLR, SVM, kNN, and RF), the stacking algorithm was able to reduce the RMSE by 25.2%, 24.9%, 22.8%, and 18.7%, respectively. To conclude, by combining ICESat-2 and continuous Sentinel-2 images, the forest canopy heights retrieved from the stacking algorithm were more consistent with the actual distribution of forest canopy heights derived from each individual algorithm, suggesting that our stacking algorithm can effectively estimate and map the spatial distribution of the forest canopy height and provide a reference for forest dynamic change monitoring and forest management.