remote sensing An Improved Generalized Hierarchical Estimation Framework with Geostatistics for Mapping Forest Parameters and Its Uncertainty: A Case Study of Forest Canopy Height

: Forest canopy height is an essential parameter in estimating forest aboveground biomass (AGB), growing stock volume (GSV), and carbon storage, and it can provide necessary information in forest management activities. Light direction and ranging (LiDAR) is widely used for estimating canopy height. Considering the high cost of acquiring LiDAR data over large areas, we took a two-stage up-scaling approach in estimating forest canopy height and aimed to develop a method for quantifying the uncertainty of the estimation result. Based on the generalized hierarchical model-based (GHMB) estimation framework, a new estimation framework named RK-GHMB that makes use of a geostatistical method (regression kriging, RK) was developed. In this framework, the wall-to-wall forest canopy height and corresponding uncertainty in map unit scale are generated. This study was carried out by integrating plot data, sampled airborne LiDAR data, and wall-to-wall Ziyuan-3 satellite (ZY3) stereo images. The result shows that RK-GHMB can obtain a similar estimation accuracy ( r = 0.92, MAE = 1.50 m) to GHMB ( r = 0.92, MAE = 1.52 m) with plot-based reference data. For LiDAR-based reference data, the accuracy of RK-GHMB ( r = 0.78, MAE = 1.75 m) is higher than that of GHMB ( r = 0.75, MAE = 1.85 m). The uncertainties for all map units range from 1.54 to 3.60 m for the RK-GHMB results. The values change between 1.84 and 3.60 m for GHMB. This study demonstrates that this two-stage up-scaling approach can be used to monitor forest canopy height. The proposed RK-GHMB approach considers the spatial autocorrelation of neighboring data in the second modeling stage and can achieve a higher accuracy.


Introduction
Forests play an essential role in climate change, ecological balance, and the carbon cycle [1][2][3]. Forest canopy height is the main structural parameter of forests. The accurate estimation of canopy height can be beneficial for the modeling of ecosystem services, forest biomass, or other forest parameters [4][5][6]. Traditional measurements by ground plots take lots of time and have a high labor cost, as well as being hard to perform, especially in dense forest areas [7,8]. Moreover, ground plot-based measurements can provide limited information about forest resources in terms of spatial coverage, because large areas cannot be surveyed due to topographic factors, the climate, or other reasons [9,10].
In the last few years, more and more studies have been conducted on estimating forest canopy height with remote sensing data [11][12][13][14]. Space-borne remote sensors can obtain data over spatially continuous large areas with a low cost. However, space-borne optical images, which are usually used to estimate the volume and aboveground biomass of forests [15][16][17], are not sensitive to forest vertical structure. Among the numerous remote

Plot Data
The field campaign was carried out in January and February 2018. A total of 54 plots (20 m × 20 m) were collected in the Jiepai and Dongsheng sub-farms ( Figure 2). During this field campaign, all trees with a minimum diameter at breast height (DBH) of 5 cm in the sample plots were measured. The tree height was measured using a handheld laser altimeter and DBH was measured with a tape. The canopy height of each plot was calculated using Equation (1).
where is the canopy height of a plot; ℎ is the height of the th individual tree in the plot; is the basal area of the th individual tree in the plot; and is the total number of trees in the plot.

Plot Data
The field campaign was carried out in January and February 2018. A total of 54 plots (20 m × 20 m) were collected in the Jiepai and Dongsheng sub-farms ( Figure 2). During this field campaign, all trees with a minimum diameter at breast height (DBH) of 5 cm in the sample plots were measured. The tree height was measured using a handheld laser altimeter and DBH was measured with a tape. The canopy height of each plot was calculated using Equation (1).
where H is the canopy height of a plot; h i is the height of the ith individual tree in the plot; g i is the basal area of the ith individual tree in the plot; and N is the total number of trees in the plot.

ALS Data Acquisition and Processing
Airborne LiDAR data were acquired to cover the whole study area in January 2018 using a CAF-LiCHy system [52] at an altitude of 1000 m above ground level. The specific parameters of the LiDAR sensor and imaging platform are shown in Table 1. The average point density was 8.51 pts/m 2 . Some noise points were deleted by setting the highest and lowest thresholds. The non-ground and ground points were classified by the Kraus filter algorithm. DTM derived from LiDAR data (DTMLiDAR) in a 1 m × 1 m resolution (grid cell) was produced by calculating the mean height of ground points within one grid cell. The non-ground point clouds were then normalized by this DTM. After that, LiDAR heightrelated variables, including the 50th/75th/95th percentiles of height named Hp50/Hp75/Hp95, were calculated based on the normalized point clouds. Here, we selected 37 plots and these LiDAR variables for modeling and estimating forest canopy height (HLi) covering the entire study area ( Figure 3). The left 17 plots were used for validation with an accuracy of = 0.95 and = 1.56 m.

ALS Data Acquisition and Processing
Airborne LiDAR data were acquired to cover the whole study area in January 2018 using a CAF-LiCHy system [52] at an altitude of 1000 m above ground level. The specific parameters of the LiDAR sensor and imaging platform are shown in Table 1. The average point density was 8.51 pts/m 2 . Some noise points were deleted by setting the highest and lowest thresholds. The non-ground and ground points were classified by the Kraus filter algorithm. DTM derived from LiDAR data (DTM LiDAR ) in a 1 m × 1 m resolution (grid cell) was produced by calculating the mean height of ground points within one grid cell. The non-ground point clouds were then normalized by this DTM. After that, LiDAR height-related variables, including the 50th/75th/95th percentiles of height named Hp50/Hp75/Hp95, were calculated based on the normalized point clouds. Here, we selected 37 plots and these LiDAR variables for modeling and estimating forest canopy height (H Li ) covering the entire study area ( Figure 3). The left 17 plots were used for validation with an accuracy of r = 0.95 and MAE = 1.56 m.

ZY3 Stereo Images and Processing
The ZY3 surveying satellite is a high-precision stereo mapping satellite of China that was launched on 9 January 2012. Two push-broom imaging sensors were carried on the satellite for acquiring multispectral and panchromatic images. The multispectral sensor can obtain a 4-band (blue, green, red and near infrared) image with a spatial resolution of 5.8 m × 5.8 m. The panchromatic sensor can obtain nadir-view images (2.1 m × 2.1 m), forward images (3.5 m × 3.5 m), and backward images (3.5 m × 3.5 m).
The ZY3 stereo images for the test site were acquired on 10 March 2018. The DSM derived from ZY3 data (DSMZY3) was generated by a dense image matching method from three panchromatic stereo images; the spatial resolution was 5 m × 5 m. Then, one CHM with a 5 m × 5 m resolution was produced by taking the difference between DSMZY3 and DTMLiDAR, which was sampled to 5 m × 5 m by a mean filter from its 1 m×1 m format. This CHM was further sampled to a grid size of 20 m × 20 m, the same size as our forest plot. The resulting CHM was named CHMZY3 and used as an independent variable for fitting the model in the second stage of RK-GHMB.

Overview
We used the canopy height derived from full-coverage LiDAR data as a baseline to see if the simulated LiDAR sampling data in strips combined with ZY3 stereo images and ground plot data could achieve a similar accuracy. We also wanted to confirm whether RK could be used to promote the accuracy and reduce the uncertainty of the estimation result. Figure 4 shows the workflow for this study.
Twelve strips of LiDAR data were sampled, assuming that only these areas were available for obtaining LiDAR data in a practical application. The selected strips that covered the 37 train plots are shown in Figure 5. Each strip was 600 m wide and the distance between the center lines of these strips was 1500 m. The simulated strips covered about 65% of the study area. A total of 37 field plots were selected as training data to build the first stage regression model combined with LiDAR data named model F. A total of 3000 pixels of canopy height estimated by model F were used to build the second stage model. The left 17 field plots and the 1000 pixels randomly selected from HLi were used as two kinds of validation data, named plot-based and LiDAR-based reference data, respectively. It is important to emphasize that these LiDAR-reference validation data are different from

ZY3 Stereo Images and Processing
The ZY3 surveying satellite is a high-precision stereo mapping satellite of China that was launched on 9 January 2012. Two push-broom imaging sensors were carried on the satellite for acquiring multispectral and panchromatic images. The multispectral sensor can obtain a 4-band (blue, green, red and near infrared) image with a spatial resolution of 5.8 m × 5.8 m. The panchromatic sensor can obtain nadir-view images (2.1 m × 2.1 m), forward images (3.5 m × 3.5 m), and backward images (3.5 m × 3.5 m).
The ZY3 stereo images for the test site were acquired on 10 March 2018. The DSM derived from ZY3 data (DSM ZY3 ) was generated by a dense image matching method from three panchromatic stereo images; the spatial resolution was 5 m × 5 m. Then, one CHM with a 5 m × 5 m resolution was produced by taking the difference between DSM ZY3 and DTM LiDAR , which was sampled to 5 m × 5 m by a mean filter from its 1 m×1 m format. This CHM was further sampled to a grid size of 20 m × 20 m, the same size as our forest plot. The resulting CHM was named CHM ZY3 and used as an independent variable for fitting the model in the second stage of RK-GHMB.

Overview
We used the canopy height derived from full-coverage LiDAR data as a baseline to see if the simulated LiDAR sampling data in strips combined with ZY3 stereo images and ground plot data could achieve a similar accuracy. We also wanted to confirm whether RK could be used to promote the accuracy and reduce the uncertainty of the estimation result. Figure 4 shows the workflow for this study.
Twelve strips of LiDAR data were sampled, assuming that only these areas were available for obtaining LiDAR data in a practical application. The selected strips that covered the 37 train plots are shown in Figure 5. Each strip was 600 m wide and the distance between the center lines of these strips was 1500 m. The simulated strips covered about 65% of the study area. A total of 37 field plots were selected as training data to build the first stage regression model combined with LiDAR data named model F. A total of 3000 pixels of canopy height estimated by model F were used to build the second stage model. The left 17 field plots and the 1000 pixels randomly selected from H Li were used as two kinds of validation data, named plot-based and LiDAR-based reference data, respectively. It is important to emphasize that these LiDAR-reference validation data are different from those 3000 pixels used to build the second stage model. For the convenience those 3000 pixels used to build the second stage model. For the convenience of subsequent descriptions, we named the datasets used in the first and second regression models as S and Sa, respectively.

Regression Model
Ordinary Least Squares (OLS) regression was used in the first stage and second stage regression models for applying GHMB to our datasets. The fitted regression models were those 3000 pixels used to build the second stage model. For the convenience of subsequent descriptions, we named the datasets used in the first and second regression models as S and Sa, respectively.

Regression Model
Ordinary Least Squares (OLS) regression was used in the first stage and second stage regression models for applying GHMB to our datasets. The fitted regression models were

GHMB Regression Model
Ordinary Least Squares (OLS) regression was used in the first stage and second stage regression models for applying GHMB to our datasets. The fitted regression models were Remote Sens. 2022, 14, 568 7 of 17 named F and G, respectively, in our study. The forms of the equation can be expressed as follows: where α 0 and β 0 are the intercepts; α j and β j are the coefficients; p is the count of variables; ε i and ξ i are the residual errors; and i is the train sample index.
In model F, Y F i is the forest canopy height of each plot selected for constructing the model, and 37 plots in total were used for modeling. X is the LiDAR variables; here, p can change from 1 to 3 according to the decisions made in the selection of the variables.
In model G, Y G i is the forest canopy height of one pixel randomly selected from the LiDAR product, and 3000 pixels in total were used for modeling. X is the CHM ZY3 , here p = 1.

Uncertainties Estimation
Saarela et al. [10] gave a method to describe the uncertainties of the two-stage estimation results in terms of the root mean square error (RMSE). The form of this method can be written as follows: There are two parts in Equation (4) that represent the uncertainty caused by the coefficients and residual error. Here,Ŷ G i is the result estimated by model G at map unit i; Z i is a (p + 1) length vector of partial derivatives of model G with respect to α Sa ; p is the count of the variables; and Z Sa is a matrix of partial derivatives of model G with respect to α Sa based on dataset Sa. V(ξ i ) is the variance of the individual random error ξ i when the heteroskedasticity of residuals exists; otherwise, V(ξ i ) is the variance of the residuals.
C is the n × n covariance matrix of residuals for the dataset Sa. The diagonal elements are estimated by V(ξ i ). In the case of autocorrelation [53], the off-diagonal elements are recalculated as follows: Cov The estimated spatial correlation ρ ij of the two residuals is calculated as follow: γ d ij is an exponential semi-variogram that was calculated using a variofit() function in the R package geoR [54].

RK-GHMB Regression Model
The regression model of the first stage of RK-GHMB is the same as GHMB, e.g., Equation (2). However, we use the RK model instead of the OLS model in the second stage of GHMB.
RK is a geostatistical prediction technique that combines the regression result and kriging result of the model residuals [55]. This technique can remove the trends in the estimates made by the regression model and has been proven to be able to mitigate the saturation problem in estimating forest attributes [56]. RK can be briefly described as the sum of the regression trend and the residuals trend [57]. The form of RK can be expressed as follows: where Y G i is the estimated result of model G, and R kriging is the residual result interpolated by ordinary kriging. R kriging can be calculated as follows: where j is the index of model G residuals, ξ j is the residual value at j index, n is the number of total residuals, and ω j is the kriging weight; the total weights must add up to one. The value of ω j can be derived as follows: where ω is a n length vector of ω j kriging weights; C is the n × n covariance matrix of residuals; and c i is an n length vector of covariance for the map unit i.

Uncertainties Estimation
Cressie et al. [58] gave the equation for estimating the uncertainty of the RK result as Equation (13). The first part is associated with the error of estimating the regression trend, and the second part presents the kriging variance of residuals. However, this method cannot satisfy the two-stage upscale processing. We modify it and obtain the equation for this process as Equation (14). Here, Cov(α Sa ) is replaced by Equation (5), C 0 + C 1 is the variance of residuals, and V(ξ i ) is equal to it when the heteroskedasticity of the residuals does not exist.

Accuracy Assessment
The ZY3 estimated canopy heights generated from GHMB and RK-GHMB were compared to the plot-based and LiDAR-based reference data with the accuracy reported using r and MAE. where Y is the group heights of the reference data;Ŷ is the group of predicted heights; y i is the height of the ith data in Y;ŷ i is the height of the ith data inŶ; and N is the number of reference data.

Forest Canopy Height Estimation Result of GHMB
Field plots and LiDAR-derived variables were used to establish the first-stage model F through the gnls() function in the R package nlme [59]. The regression model F is shown as Equation (17). Hp95, with the highest correlation (p = 0.87), was selected as the variable for model F.
In these sampled LiDAR strips, 3000 pixels of LiDAR-derived canopy height were used to build the second-stage regression model G. The equation of model G is shown as Equation (18). Then, the wall-to-wall CHM ZY3 raster layer with a high correlation (p = 0.80) was applied in Equation (18) to generate a ZY3-estimated canopy height map. To describe the heteroskedasticity, we fitted a nonlinear power model to estimate the individual error variance for each ZY3-estimated value. The equation fitted is shown as Equation (19) in Table 2.

Forest Canopy Height Estimation Result of RK-GHMB
A variogram of residuals generated by model G was calculated. An exponential model was selected to fit this variogram ( Figure 6). The basic parameters of the fitted exponential model for ordinary kriging are shown in Table 3. The ratio of the nugget and sill in the fitted model (C 0 /(C 0 + C 1 )) represents the spatial dependence structure. The higher the ratio is, the more variations are determined by random effects [60]. Usually, a ratio under 25% means that there is a strong spatial dependence structure. If the ratio is between 25% and 75%, the spatial dependence structure is moderate, while it is weak if the ratio is higher than 75%. In our study, the ratio is 33.2%, which means that the residuals had a moderate spatial dependence structure. The practical range is 147.27 m, which means that the RK model could have a limited effect out of this distance. Based on Figure 6 and Table 3, a residual-kriging map was generated. The distribution of canopy height (Figure 7b) based on the RK-GHMB model was acquired by combining the residual-kriging map and the model G-estimated map (Figure 7a).

Forest Canopy Height Estimation Accuracy and Uncertainty Assessment
We used plots and the LiDAR-derived forest canopy height as testing reference data to compare the performance of GHMB and RK-GHMB for forest canopy height estimation. Figure 8 shows the scatter plots, and Table 4 presents the accuracy of these estimation results. The plot-based and for GHMB are 0.92 and 1.52 m, while for RK-GHMB, the values are 0.92 and 1.50 m. The LiDAR-based and for GHMB are 0.75 and 1.85 m, while for RK-GHMB, the values are 0.78 and 1.75 m. The accuracy of the two models is similar to that of the plot-based reference data. However, the RK-GHMB model has a higher and a lower compared to the LiDAR-based reference data. In order to further prove the effect of improving the accuracy, we randomly selected 1000 LiDARbased reference data 100 times in order to check the accuracy and obtain two groups of r and MAE. Then, we checked whether there was a significant difference between the two groups. Figure 9 shows the corresponding r and MAE for the two estimation frameworks. According to the Wilcoxon test, a significant difference (p < 0.05) was seen between the

Forest Canopy Height Estimation Accuracy and Uncertainty Assessment
We used plots and the LiDAR-derived forest canopy height as testing reference data to compare the performance of GHMB and RK-GHMB for forest canopy height estimation. Figure 8 shows the scatter plots, and Table 4  The accuracy of the two models is similar to that of the plot-based reference data. However, the RK-GHMB model has a higher and a lower compared to the LiDAR-based reference data. In order to further prove the effect of improving the accuracy, we randomly selected 1000 LiDARbased reference data 100 times in order to check the accuracy and obtain two groups of r and MAE. Then, we checked whether there was a significant difference between the two groups. Figure 9 shows the corresponding r and MAE for the two estimation frameworks. According to the Wilcoxon test, a significant difference (p < 0.05) was seen between the

Forest Canopy Height Estimation Accuracy and Uncertainty Assessment
We used plots and the LiDAR-derived forest canopy height as testing reference data to compare the performance of GHMB and RK-GHMB for forest canopy height estimation. Figure 8 shows the scatter plots, and Table 4  The accuracy of the two models is similar to that of the plot-based reference data. However, the RK-GHMB model has a higher r and a lower MAE compared to the LiDAR-based reference data. In order to further prove the effect of improving the accuracy, we randomly selected 1000 LiDAR-based reference data 100 times in order to check the accuracy and obtain two groups of r and MAE. Then, we checked whether there was a significant difference between the two groups. Figure 9 shows the corresponding r and MAE for the two estimation frameworks. According to the Wilcoxon test, a significant difference (p < 0.05) was seen between the accuracy of GHMB and RK-GHMB. This comparison shows that the additional prediction of residuals by kriging can increase the estimation accuracy.     accuracy of GHMB and RK-GHMB. This comparison shows that the additional prediction of residuals by kriging can increase the estimation accuracy.    The uncertainty maps of GHMB and RK-GHMB are shown as Figure 10. The uncertainties estimated by GHMB ranged from 1.84 to 3.60 m. However, the uncertainties estimated using RK-GHMB ranged from 1.54 to 3.60 m. In Figure 10b, the RMSE of most map units seems to be lower than that in Figure 10a; this means that RK-GHMB can obtain a higher accuracy in forest height estimation with a lower uncertainty. From Figure 10b, we can see that the pixels with a higher uncertainty reduction are located near the twelve sampled LiDAR data strips.
The uncertainty maps of GHMB and RK-GHMB are shown as Figure 10. The uncertainties estimated by GHMB ranged from 1.84 to 3.60 m. However, the uncertainties estimated using RK-GHMB ranged from 1.54 to 3.60 m. In Figure 10b, the RMSE of most map units seems to be lower than that in Figure 10a; this means that RK-GHMB can obtain a higher accuracy in forest height estimation with a lower uncertainty. From Figure 10b, we can see that the pixels with a higher uncertainty reduction are located near the twelve sampled LiDAR data strips.
The uncertainties for the selected 1000 LiDAR reference data are shown in Figure 11. From Figure 11, we can clearly see that RK-GHMB can obtain a lower RMSE value, and the uncertainties have a similar changing trend, with the RMSE being larger at both higher and lower predicted values.

Discussion
In this work, we applied a two-stage up-scaling approach to produce a forest canopy height map. As in many previous studies, LiDAR data were used as a bridge, combining The uncertainties for the selected 1000 LiDAR reference data are shown in Figure 11. From Figure 11, we can clearly see that RK-GHMB can obtain a lower RMSE value, and the uncertainties have a similar changing trend, with the RMSE being larger at both higher and lower predicted values. The uncertainty maps of GHMB and RK-GHMB are shown as Figure 10. The uncertainties estimated by GHMB ranged from 1.84 to 3.60 m. However, the uncertainties estimated using RK-GHMB ranged from 1.54 to 3.60 m. In Figure 10b, the RMSE of most map units seems to be lower than that in Figure 10a; this means that RK-GHMB can obtain a higher accuracy in forest height estimation with a lower uncertainty. From Figure 10b, we can see that the pixels with a higher uncertainty reduction are located near the twelve sampled LiDAR data strips.
The uncertainties for the selected 1000 LiDAR reference data are shown in Figure 11. From Figure 11, we can clearly see that RK-GHMB can obtain a lower RMSE value, and the uncertainties have a similar changing trend, with the RMSE being larger at both higher and lower predicted values.

Discussion
In this work, we applied a two-stage up-scaling approach to produce a forest canopy height map. As in many previous studies, LiDAR data were used as a bridge, combining

Discussion
In this work, we applied a two-stage up-scaling approach to produce a forest canopy height map. As in many previous studies, LiDAR data were used as a bridge, combining plot data and optical remote sensing data. Considering the high cost of acquiring and processing LiDAR data, this combination can reduce the use of LiDAR data, effectively reducing the cost of estimating the forest parameters of large areas. Although the LiDARestimated forest parameters still contain some errors, more and more studies [19,24,28,29,61] prefer to use these products for plot data in model building and validation. As previous studies have claimed [28,62], this operation can increase the sample size and distribution range of sample data.
In this study, we used a geostatistical model named RK in the second stage of GHMB. The results of model RK has combined the residual data interpolated by ordinary kriging and the estimated results of model G, just as Fayad et al. [19] described for the random forest model. The wall-to-wall results of model G do not take into account the spatial correlation between the selected samples, assuming it to be spatially independent. However, some of the unexplained variance in the results of model G could be due to the spatial correlation of the selected samples. This procedure of combing the kriging residuals can remove the trends caused by model G. The advantageous performance of RK were also demonstrated in the estimation of forest canopy heights [19,29], AGB [63][64][65], and other forest parameters [66]. In previous studies, the range was almost 200 m and 4421-4823 m in the studies of Hudak et al. [29] and Fayad et al. [19], respectively. Mauro et al. [67] revealed that the spatial correlation range of residuals is reduced the most when independent variables are highly correlated with dependent variables. In our study, the CHM ZY3 raster layer has a high correlation (p = 0.80) with the LiDAR-estimated forest canopy height. This can be used to explain why the range of residuals derived by model G is so small (147.27 m). This value indicates that the kriging residual data have little influence on the RK result and limit the performance of promoting the accuracy out of this distance. The distance of the sampled LiDAR strips was 900 m, which was larger than the spatial auto-correlation range. This can explain why there is only a slight improvement in RK-GHMB compared to GHMB. There are also some limitations in the data. For example, there is a time interval of two months between plot data and ZY3 stereo images. There is a study [68] showing that the tree height grows 3.6 m a year when it is 9.6 m, and 0.8 m a year for 13.3 m by taking Eucalyptus as an example. According to the LiDAR-estimated result in our study area, the forest canopy height was almost larger than 9 m and around 15 m, which means that the tree height does not grow very much in a short period of time. However, the simultaneous acquisition of the image data should yield a higher estimation accuracy. The data size of the sample plot is relatively small, and the measured sample plots are also concentrated on Eucalyptus and Chinese Fir species. The limited tree species may affect the final estimation results. Especially when estimating biomass or other forest parameters.
In this study, we estimated the forest canopy height and the corresponding uncertainty. Uncertainties were accounted for the models used in the up-scaling processing. GHMB was used to estimate the uncertainties caused by model F and model G, while RK-GHMB was used to estimate the uncertainties caused by model F and model RK. As shown in Figure 11, the RMSE was larger at both higher and lower predicted values, whether GHMB or RK-GHMB were used for estimating the uncertainties. This phenomenon was also seen by Gu et al. [61], who claimed that shorter and higher trees had higher errors. In addition, the uncertainties of the map units near the sampled LiDAR strips were much smaller in the RK-GHMB-estimated result. This phenomenon was also seen by Fayad et al. [19], who claimed that the uncertainty of the canopy height seems to be correlated with the location of the reference dataset. Both the models and validation results tell us that the accuracy increment from RK-GHMB to GHMB is determined by the degree of spatial autocorrelation between the map units of the first-stage estimation results. The more geocorrelation there is, the more improvement there will be. If there is no spatial correlation, the RK-GHMB should shrink to GHMB. The uncertainty of the final result is multifaceted, affected by the sampling error and measurement error caused by differences in measuring instruments or observers. However, only the uncertainty caused by the model coefficients and residuals can be obtained under this estimation framework. Different from using uncertainty based on pixels, obtaining the forest parameter results and their uncertainty for small-area (stand-level) estimation will be our next research direction.
As a general estimation framework, the proposed RK-GHMB can also be used to map AGB, GSV, and other forest parameters together with their uncertainties in the scale of mapping units (the same size as a ground forest plot) on the condition that forest plot data (sampling data), airborne data (sampling data), and space or airborne remote sensing data (full coverage data) exist. In order to achieve better performance with the proposed RK-GHMB, firstly, the remote sensing features derived from airborne data should have good correlation with the target parameters computed from ground plot data; this can be achieved by collecting high-density LiDAR sampling data using UAV or a humancontrolled airplane. If the mapping area has archived airborne LiDAR data and thus a high-resolution DTM product, it can be a good solution to acquire low-cost airborne stereo images; the DSM and relevant features derived from these will normally have a good correlation with the forest parameters. Secondly, the remote sensing data fully covering the target mapping area should be much cheaper to acquire than the remote sensing data used in the first stage, while the wall-to-wall remote sensing features used in the second stage should be less sensitive to the target forest parameters. Here, space-borne stereo images, low-density air-borne LiDAR, or stereo images are some suitable data sources.

Conclusions
We proposed a new forest parameter and uncertainty estimation framework, named RK-GHMB, by considering the spatial autocorrelation of neighboring data in the second modeling stage using the RK model. We validated its performance, taking forest canopy height mapping as one application case. The results show that the RK-GHMB can achieve a higher forest canopy height estimation accuracy than GHMB and can model the estimation uncertainty for all map units with an improved performance when taking GHMB as a baseline model. We successfully demonstrated that the RK-GHMB model is a useful solution for mapping forest parameters and their uncertainty in large areas by integrating ground plot data, sampled LiDAR data, and wall-to-wall ZY3 stereo images. Data Availability Statement: The field data are not publicly available due to privacy of the private landowners.

Conflicts of Interest:
The authors declare no conflict of interest.