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

Abstract

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.

1. 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 sensing technologies available, light direction and ranging (LiDAR) technology has a strong advantage in canopy height extraction due to its ability to penetrate and obtain the vertical structure information of forests. In particular, the airborne laser scanning (ALS) system, has been widely used in canopy height estimation research [18,19,20]. However, the high cost of data acquisition restricts its application in a wide range of spatial distributions [21,22,23]. The most effective and lowest-cost way to map forest parameters over large areas may be to combine ground plots, ALS data obtained by sampling, and wall-to-wall space-borne remote sensing data [24]. Some studies have shown that the two-stage up-scaling approach is practical and acceptable in mapping forest parameters [25,26,27,28]. Regression kriging (RK), as a geostatistical method, is often used in this approach to improve the accuracy of estimation [19,29,30]. Vegetation indices, reflectance, or texture parameters extracted from space-borne optical images are common features used to estimate forest parameters [31,32,33,34,35]. Different from these features, the stereo images acquired from space-borne remote sensors can generate digital surface models (DSM). Canopy height model (CHM) and other features generated from DSM have been commonly used for mapping canopy height [36,37,38,39].

In order to evaluate the accuracy of the mapping results, commonly used uncertainty metrics are global correlation and root mean square error between the limited field reference data and the estimated value. However, the accuracy often varies in spatial distribution depending on the complexity of landscape, density of sampled data, accuracy of the remote sensed data used, etc. [40,41]. Pixel-based uncertainty estimation can reflect the accuracy of estimation results in spatial distribution. This is very useful and important for users to judge whether these estimation results are satisfied or fitted enough in further decision making. There have been many studies focused on the uncertainty of estimation results in forest parameter estimation [19,28,42,43,44,45,46,47]. For example, García et al. [20] combined plot data, LiDAR, synthetic aperture radar and Landsat data for estimating forest canopy height and its uncertainty by a bootstrapping method. Gonzalez et al. [48] estimated forest carbon and its uncertainty by Monte Carlo analysis. Lang et al. [49] estimated the forest canopy height by a probabilistic deep learning method. Among these research studies, few have considered the uncertainty in the two-stage up-scaling process [44,46,48], especially for forest canopy height [20,49]. Besides these, a method named hierarchical model-based (HMB) estimation was proposed by Saarela et al. [50]. Saarela et al. [50] pointed out that neglecting the uncertainty associated with one of the two models involved might lead to the severe underestimation of the uncertainty. HMB assumes that the estimation models involved must have a homogeneous residual variance and restricted linear models. Saarela et al. [51] further developed a new method named the generalized hierarchical model-based (GHMB) estimation framework based on HMB. GHMB can overcome some shortcomings of HMB and can be used when the models involved have non-homogeneous variance. Saarela et al. [10] improved the GHMB method further so that it can use nonlinear models.

Although the RK model has better performance than some regression models when spatial autocorrelation exists between units for parameter estimation, it has not been used in the second stage of either HMB or GHMB. Assuming that we can improve the performance of GHMB by replacing the regression model with the RK model in the second estimation stage, one general forest parameter estimation and uncertainty estimation framework named RK-GHMB has been developed in this paper. The validation of RK-GHMB was carried out through forest canopy height estimation by the combined use of ground plot data, sampled LiDAR data, and Ziyuan-3 satellite (ZY3) stereo images.

2. Materials and Methods

2.1. Study Area

Gaofeng forest farm is situated in Nanning, Guangxi Zhuang Autonomous Region, China. This region has a humid subtropical monsoon climate with an average annual temperature of 21 °C and an average annual precipitation of 1200–1500 mm. Typical tree species in this farm are Eucalyptus, Chinese Fir, Masson Pine, and other broad-leaved tree species. The study area is located in two sub-farms named Jiepai and Dongsheng within the Gaofeng forest farm with an area of 52 km² (Figure 1). The slope ranges from 0° to 69.7°, and the elevation ranges from 77 to 463 m.

2.2. Data

2.2.1. 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).

$$H=\frac{{{\displaystyle \sum }}_{i=1}^N{h}_i{g}_i}{{{\displaystyle \sum }}_{i=1}^N{g}_i}$$

(1)

where $H$ is the canopy height of a plot; $h_i$ is the height of the $i\mathrm{th}$ individual tree in the plot; $g_i$ is the basal area of the $i\mathrm{th}$ individual tree in the plot; and $N$ is the total number of trees in the plot.

2.2.2. 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. Parameters of the LiDAR sensor and imaging platform.

Parameters	Value	Parameters	Value
Platform	Tecnam P2006T	Flying height (m)	1000
Laser beam divergence (m·rad)	0.5	Speed (Km·h−1)	180
Laser wavelength (nm)	1550	Vertical accuracy (cm)	15
Scan angle (°)	±30	Average point density (points·m−2)	8.51
Laser pulse repetition rate (kHz)	400	Pulse length (ns)	3

. The average point density was 8.51 pts/m². 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.

2.2.3. 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.

2.3. Methods

2.3.1. 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 of subsequent descriptions, we named the datasets used in the first and second regression models as S and Sa, respectively.

2.3.2. 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 named F and G, respectively, in our study. The forms of the equation can be expressed as follows:

$$Y_{F_i}={\beta }_0+{\displaystyle \sum }_{j=1}^p{\beta }_j{X}_i+{\epsilon }_i$$

(2)

$$Y_{G_i}={\alpha }_0+{\displaystyle \sum }_{j=1}^p{\alpha }_j{X}_i+{\xi }_i$$

(3)

where ${\alpha }_0$ and ${\beta }_0$ are the intercepts; ${\alpha }_j$ and ${\beta }_j$ are the coefficients; $p$ is the count of variables; ${\epsilon }_i$ and ${\xi }_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:

$$RMSE\left({\widehat{Y}}_{G_i}\right)=\sqrt{Z_i^{T}Cov\left({\widehat{\alpha }}_{Sa}\right)Z_i+V\left({\xi }_i\right)}$$

(4)

$$Cov\left({\widehat{\alpha }}_{Sa}\right)={\left(Z_{Sa}^T·C^{-1}·Z_{Sa}\right)}^{-1}+{\left(Z_{Sa}^T·C^{-1}·Z_{Sa}\right)}^{-1}·Z_{Sa}^T·C^{-1}·Cov\left({\widehat{Y}}_{F_{Sa}}\right)·C^{-1}·Z_{Sa}·{\left(Z_{Sa}^T·C^{-1}·Z_{Sa}\right)}^{-1}$$

(5)

There are two parts in Equation (4) that represent the uncertainty caused by the coefficients and residual error. Here, ${\widehat{Y}}_{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 ${\alpha }_{Sa}$; p is the count of the variables; and $Z_{Sa}$ is a matrix of partial derivatives of model G with respect to ${\alpha }_{Sa}$ based on dataset Sa. $V\left({\xi }_i\right)$ is the variance of the individual random error ${\xi }_i$ when the heteroskedasticity of residuals exists; otherwise, $V\left({\xi }_i\right)$ is the variance of the residuals.

$C$is the $n\times n$covariance matrix of residuals for the dataset Sa. The diagonal elements are estimated by $V\left({\xi }_i\right)$. In the case of autocorrelation [53], the off-diagonal elements are recalculated as follows:

$$Cov\left({\xi }_i,{\xi }_j\right)=\sqrt{V\left({\xi }_i\right)·V\left({\xi }_j\right)}·{\rho }_{ij}$$

(6)

The estimated spatial correlation ${\rho }_{ij}$ of the two residuals is calculated as follow:

$${\rho }_{ij}=1-\frac{\gamma \left(d_{ij}\right)}{C_0+C_1}$$

(7)

$\gamma \left(d_{ij}\right)$ is an exponential semi-variogram that was calculated using a variofit() function in the R package geoR [54].

2.3.3. 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:

$$Y_{R{K}_i}=Y_{G_i}+R_{kriging}$$

(8)

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:

$$R_{kriging}={\displaystyle \sum }_{j=1}^n{\omega }_j·{\xi }_j$$

(9)

$${\displaystyle \sum }_{j=1}^n{\omega }_j=1$$

(10)

where $j$ is the index of model G residuals, ${\xi }_j$ is the residual value at $j$ index, $n$ is the number of total residuals, and ${\omega }_j$ is the kriging weight; the total weights must add up to one. The value of ${\omega }_j$ can be derived as follows:

$$\omega =C^{-1}·c_i$$

(11)

$$c_i=\left(\begin{array}{c}\begin{array}{c}Cov\left({\xi }_i,{\xi }_1\right)\\ Cov\left({\xi }_i,{\xi }_2\right)\end{array}\\ ⋮\\ \begin{array}{c}Cov\left({\xi }_i,{\xi }_{n-1}\right)\\ Cov\left({\xi }_i,{\xi }_n\right)\end{array}\end{array}\right)$$

(12)

where $\omega$ is a $n$ length vector of ${\omega }_j$ kriging weights; $C$ is the$n\times 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\left({\widehat{\alpha }}_{Sa}\right)$ is replaced by Equation (5), $C_0+C_1$ is the variance of residuals, and $V\left({\xi }_i\right)$ is equal to it when the heteroskedasticity of the residuals does not exist.

$${\sigma }_{{\widehat{Y}}_{R{K}_i}}^2={\left(Z_i-Z_{Sa}^T·C^{-1}·c_i\right)}^T·Cov\left({\widehat{\alpha }}_{Sa}\right)·\left(Z_i-Z_{Sa}^T·C^{-1}·c_i\right)+\left(C_0+C_1\right)-c_i^T·C^{-1}·c_i$$

(13)

$$RMSE({\stackrel{⌢}{Y}}_{R{K}_i})=\sqrt{{\sigma }_{{\stackrel{⌢}{Y}}_{R{K}_i}}^2}=\sqrt{\begin{array}{l}(Z_i-Z_{Sa}^T·C^{-1}·c_i{)}^T\\ ·\left\{{\left(Z_{Sa}^T·C^{-1}·Z_{Sa}\right)}^{-1}+{\left(Z_{Sa}^T·C^{-1}·Z_{Sa}\right)}^{-1}·Z_{Sa}^T·C^{-1}·Cov({\stackrel{⌢}{Y}}_{F_{Sa}})·C^{-1}·Z_{Sa}·{\left(Z_{Sa}^T·C^{-1}·Z_{Sa}\right)}^{-1}\right\}\\ ·\left(Z_i-Z_{Sa}^T·C^{-1}·c_i\right)+\left\{V({\xi }_i)-c_i^T·C^{-1}·c_i\right\}\end{array}}$$

(14)

2.3.4. 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$.

$$r=\frac{Cov\left(Y,\widehat{Y}\right)}{\sqrt{Var\left(Y\right)·Var\left(\widehat{Y}\right)}}$$

(15)

$$MAE=\frac{1}{N}{\displaystyle \sum }_{i=1}^N\lfloor y_i-{\widehat{y}}_i\rfloor$$

(16)

where $Y$ is the group heights of the reference data; $\widehat{Y}$ is the group of predicted heights; $y_i$ is the height of the $i\mathrm{th}$ data in $Y$; ${\widehat{y}}_i$ is the height of the $i\mathrm{th}$ data in $\widehat{Y}$; and $N$ is the number of reference data.

3. Results

3.1. 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. Fitted forms for model F, model G and the variance for $\xi$.

Model Name	Model Forms		R2	RMSE
F	$$Y_F=2.1315+0.7706\times Hp95$$	(17)	0.75	1.81
G	$$Y_G=5.7147+0.6007\times CH{M}_{ZY3}$$	(18)	0.64	2.38
$$V\left({\xi }_i\right)$$	$$V\left({\xi }_i\right)=24.9245-2.2393\times Y_{G_i}+0.0577\times Y_{G_i}^2$$	(19)	0.81	1.01

.

3.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. Parameters of the fitted exponential model for model G residuals.

Residuals Source	Model	Nugget (C0)	Partial Sill (C1)	Range (m)	Ratio (%)
Model G	exponential	1.73	3.48	147.27	33.2

. The ratio of the nugget and sill in the fitted model (C₀/(C₀ + C₁)) 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).

3.3. 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. Models were evaluated by plot- and LiDAR-based reference data.

	Plot-Based Reference		LiDAR-Based Reference
Models	$r$	$MAE$	$r$	$MAE$
GHMB	0.92	1.52	0.75	1.85
RK-GHMB	0.92	1.50	0.78	1.75

presents the accuracy of these estimation results. The plot-based $r$ and $MAE$ 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 $r$ and $MAE$ 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 $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.

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.

4. 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 LiDAR-estimated 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 geo-correlation 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 human-controlled 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.

5. 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.