Use of Remotely Sensed Data to Enhance Estimation of Aboveground Biomass for the Dry Afromontane Forest in South-Central Ethiopia

Periodic assessment of forest aboveground biomass (AGB) is essential to regulate the impacts of the changing climate. However, AGB estimation using field-based sample survey (FBSS) has limited precision due to cost and accessibility constraints. Fortunately, remote sensing technologies assist to improve AGB estimation precisions. Thus, this study assessed the role of remotely sensed (RS) data in improving the precision of AGB estimation in an Afromontane forest in south-central Ethiopia. The research objectives were to identify RS variables that are useful for estimating AGB and evaluate the extent of improvement in the precision of the remote sensing-assisted AGB estimates beyond the precision of a pure FBSS. Reference AGB data for model calibration and estimation were collected from 111 systematically distributed circular sample plots (SPs) of 1000 m2 area. Independent variables were derived from Landsat-8, Sentinel-2 and PlanetScope images acquired in January 2019. The area-weighted mean and standard deviation of the spectral reflectance, spectral index and texture (only for PlanetScope) variables were extracted for each SP. A maximum of two independent variables from each image type was fitted to a generalized linear model for AGB estimation using model-assisted estimators. The results of this study revealed that the Landsat-8 model with the predictor variable of shortwave infrared band reflectance and the PlanetScope model with the predictor variable of green band reflectance had estimation efficiency of 1.40 and 1.37, respectively. Similarly, the Sentinel-2 model, which had predictor variables of shortwave infrared reflectance and standard deviation of green leaf index, improved AGB estimation with the relative efficiency of 1.68. Utilizing freely available Sentinel-2 data seems to enhance the AGB estimation efficiency and reduce cost and extensive fieldwork in inaccessible areas.


Introduction
Forests are paramount in regulating the global environment, mainly through sequestering carbon [1]. They are particularly important these days to combat the changing climate, which affects people's lives in many aspects. Due to the multiple significance of forest resources, information about the resource base, its spatial distribution and spatio-temporal changes have become a global

Description of the Study Area
The study was conducted in the Degaga-Gambo forest in south-central Ethiopia. It belongs to a state-owned enterprise, Oromia Forest and Wildlife. The study area is located on the eastern escarpment of the central rift valley of Ethiopia, in the Horn of Africa (Figure 1). It extends geographically from 38 • 45 to 38 • 56 E longitude and from 7 • 13 to 7 • 33 N latitude. The forest has an area of 14,176 ha. The altitude of the study area ranges from 2100 to 2730 m above sea level. The study area has a bimodal rainfall distribution. The main rainy season is from July to September while the short rainy season is from March to May [41]. The mean annual precipitation and temperature in the area are 1245 mm and 14.9 • C, respectively. geographically from 38°45' to 38°56' E longitude and from 7°13' to 7°33' N latitude. The forest has an area of 14,176 ha. The altitude of the study area ranges from 2100 to 2730 m above sea level. The study area has a bimodal rainfall distribution. The main rainy season is from July to September while the short rainy season is from March to May [41]. The mean annual precipitation and temperature in the area are 1245 mm and 14.9 °C, respectively.
The forest area has both natural and plantation forest types. The major species of plantation forest compartments, which are mostly found in the lower elevations, are Cupressus lucitanica, Pinus patula, Grevillea robusta and different Eucalyptus species. The natural forest has high tree species diversity. The dominant tree species observed in the natural forest include Syzygium guineense, Afrocarpus falcatus, Juniperus procera, Pitosporum viridiflorum, Maesa lanceolate, Millettia ferruginea, Croton macrostachyus and Maytenus arbutifolia. The objectives of the enterprise are the production of lumber and poles from the plantations and conserving the natural forests. The natural forests are home to a wide range of wildlife species and are sources of water for the downstream areas. Nevertheless, the forests are under severe pressure. Illegal cutting of trees and land-use change for settlement and farmland expansion are the common problems in the area.
The forest has complex vertical and horizontal structures. Besides the species diversity, there is large variability in tree height and wood basic density of the study forest. The mean (and range) of observed tree height was 13.90 m (4.90-40.10 m); while the mean (and range) of wood basic density (g cm −3 ) for tree species in the forest was 0.59 (0.43-0.98) [42].

Field Data Collection
The sampling frame was defined to include the Degaga-Gambo forest territory, which contains both the natural and plantation forest types. Circular sample plots (SPs) of 17.85 m radius aligned in a systematic grid at an interval of 1.18 km were used for field data collection ( Figure 1). One hundred and eleven plots (from the natural forests, plantation forests and other categories like clear-cut, The forest area has both natural and plantation forest types. The major species of plantation forest compartments, which are mostly found in the lower elevations, are Cupressus lucitanica, Pinus patula, Grevillea robusta and different Eucalyptus species. The natural forest has high tree species diversity. The dominant tree species observed in the natural forest include Syzygium guineense, Afrocarpus falcatus, Juniperus procera, Pitosporum viridiflorum, Maesa lanceolate, Millettia ferruginea, Croton macrostachyus and Maytenus arbutifolia. The objectives of the enterprise are the production of lumber and poles from the plantations and conserving the natural forests. The natural forests are home to a wide range of wildlife species and are sources of water for the downstream areas. Nevertheless, the forests are under severe pressure. Illegal cutting of trees and land-use change for settlement and farmland expansion are the common problems in the area. The forest has complex vertical and horizontal structures. Besides the species diversity, there is large variability in tree height and wood basic density of the study forest. The mean (and range) of Remote Sens. 2020, 12, 3335 5 of 23 observed tree height was 13.90 m (4.90-40.10 m); while the mean (and range) of wood basic density (g cm −3 ) for tree species in the forest was 0.59 (0.43-0.98) [42].

Field Data Collection
The sampling frame was defined to include the Degaga-Gambo forest territory, which contains both the natural and plantation forest types. Circular sample plots (SPs) of 17.85 m radius aligned in a systematic grid at an interval of 1.18 km were used for field data collection ( Figure 1). One hundred and eleven plots (from the natural forests, plantation forests and other categories like clear-cut, cropland, settlement and grassland cover types) were sampled from February 2018 to January 2019. Handheld global positioning system (GPS) receiver was used to navigate to the pre-defined locations of the SPs. Then, the precise coordinates of the plot centers were determined using differential GPS and global navigation satellite system (GLONASS) measurements. Two Topcon legacy-E + 40 dual-frequency receivers were used for this purpose [43]; one serving as a base station and the other as a rover field unit. The receivers record pseudo-range and carrier phase of GPS and GLONASS.
The base station was set up at Wondo Genet College of Forestry and Natural Resources campus. The Euclidean distance between the base station and the plot centers ranged between 21.70 and 57.20 km with an average distance of 41.80 km. To determine the position of the base station using precise point positioning, the GPS and GLONASS data were recorded continuously for 24 h [44]. At the plot centers, the rover was mounted on a 2.98 m carbon rod and recorded for 41.50 min on average using a one-second logging rate. The recordings were post-processed using the Magnet tools software [45]. The standard error of the post-processed planimetric plot coordinates ranged from 0.02 to 1.11 m with a mean of 0.23 m.
In each of the SPs, we recorded species names and measured DBH, i.e., the diameter of trees at 1.3 m above the ground, for all the trees with DBH ≥ 5 cm. Caliper or diameter tape was used for DBH measurement depending on tree size. Tree height measurements were carried out for 10 trees selected systematically in each of the plots using a Haglöf vertex laser 5 instrument [46]. Heights of the trees for which height was not measured were predicted using height-diameter models developed based on the sample trees [16,19,47].

Plot-Level AGB Estimation
Plot-level AGB was estimated by aggregating the predicted individual tree AGB in the respective plots. For predicting tree AGB in the natural forests, the allometric model constructed by [42] was used. This model has DBH, height and wood basic density as predictor variables. Wood basic density values were obtained from [48]. For plantation forests, tree AGB was estimated using species-specific allometric models. Accordingly, for Cupressus lusitanica, we used the model by [49] with DBH and height as predictor variables. For Eucalyptus species and Grevillea robusta, models by [50,51] were used, respectively, having DBH and height as predictor variables. The plot-level AGB data in units of kg m −2 were converted to Mg ha −1 (megagrams per hectare) since the data were collected from large plots (1000 m 2 ). The plot-level AGB values ranged from 0 to 845.70 Mg ha −1 with a mean and standard deviation of 184.35 Mg ha −1 and 155.10 Mg ha −1 , respectively.

Satellite Image Acquisition
Satellite images acquired in January 2019 were considered since this is the dry season when most of the undergrowth vegetation dries up and is easier to distinguish from the trees. This time window was also within the field inventory period. Additionally, selected images were restricted to those with cloud cover < 5%. A detail description of the images used in this study is given in Table 1.
Single tiles of each of the L8 and S2 products were downloaded from the USGS Earth Explorer website [52]. Both images were Level-1C products, which means that the images were corrected for any possible topographic and geometric errors. The processing level of the L8 image used in this study was L1-TP, which is a Level-1 precision and terrain corrected product. Besides terrain and topographic Remote Sens. 2020, 12, 3335 6 of 23 correction, radiometric correction has already been done for S2 products before delivery. The SBs used in this study (i.e., blue (B), green (G), red (R), near-infrared (NIR), shortwave infrared-1 (SWIR1) (for both L8 and S2), red-edge (RE) (only for S2)) have spatial resolutions of 30 m for L8 and 10 or 20 m for S2 (see Table 1 for details of the resolutions of individual bands).
We downloaded the PS Ortho Scene Product (Level-3B) from the Planet Explorer website [53]. Six scenes of orthorectified scaled Top of Atmosphere Radiance (at sensor) images were downloaded to cover the study area. These images contain information about the B, G, R and NIR SBs.

Image Processing and Independent Variable Definition
In the current study, we first evaluated a great number of potential candidate variables that could be useful for AGB modelling. A series of image processing techniques were applied to the satellite images to get the independent variables. First, atmospheric correction was done using the QGIS software version 3.1.0 [54] and python codes. For L8 and S2 images, the semi-automatic classification plugin (SCP) of QGIS was used for running the dark-object subtraction (DOS-1) algorithm, which removes the dark pixels that result from atmospheric scattering. The satellite images were transformed from spectral radiance to top of atmosphere reflectance values based on the conversion factors in the metadata file that comes along with the image files. However, the PS images were processed using the empirical line correction for conversion of radiance to reflectance values indicated in Equation (1): The radiances of the input images were converted to reflectance values and atmospheric correction applied since variables from multiple images were compared. In addition to variation in the sensors, the three sets of images were acquired on different dates although within 13 days of maximum gap among them. Furthermore, six scenes of the PS imagery covered the area of interest. After atmospheric correction, all the images became Level-2A products, which have pixels with surface reflectance values suitable for calculating SIs and texture variables used in this study. Atmospherically corrected SBs, which were used for creating SIs and texture variables shown in Table 2 and Table 3, respectively, were selected for this study.
It measures the degree of randomness of pixel values in an image. Entropy is inversely related to uniformity.

GLCM mean
Mean of GL distribution of the image.
GLCM variance Correlation indicates the linear dependency of GL on their neighboring pixels. d IDM and ASM stand for inverse difference moment and angular second moment, respectively. e Where p i, j is the probability of finding the GLCM relationship at cell (i, j) and is calculated as such that N−1 i, j=0 p i, j = 1, N = Number of grey levels in the image as specified by the number of levels in the quantization, V i, j = grey level value in a cell (i, j) of the image window. Table 2 shows the expressions used to derive spectral index values from each satellite image type used in this study and references to scientific evidences on the use of the indices in general and for biomass estimation in particular.
Descriptions of the GLCM image texture data derived from the PS images are presented in Table 3. Texture information of the L8 and S2 images were not used due to the coarse spatial resolutions. Sentinel Application Platform (SNAP) software version 7.0.0 [69] was used for calculating the texture variables. Processing parameters of window size of 11 × 11 pixels, angle in all directions, Remote Sens. 2020, 12, 3335 8 of 23 probability quantization with level of 128 were set to obtain the texture data used in the current study. This processing window size was set to provide an equivalent area to the field SPs.
Area-weighted mean and standard deviation (hereafter referred to as mean and standard deviation, respectively) of all the variables were extracted to each SP using QGIS. These were used as independent variables of the models constructed from each RS data type, the details of which are explained in the following sections.

Variable Selection and Model Fitting
The purpose of the AGB regression modelling was to construct models with variables from the RS data as predictors and which could be used to enhance the precision of the overall AGB estimates for the study area. For the AGB estimation, we used a model-assisted approach to inference (see details in Section 2.8) because that would allow a direct comparison of the uncertainty of the AGB estimate with similar uncertainty estimates obtained for the pure field-based estimate. In model-assisted estimation, the model form and the predictors selected for the model should be determined independent of the sample at hand [70]. In model-assisted inference, no claim of a true model is necessary. A poor choice of model form and predictors would have negative consequences in terms of efficiency [71] (p. 238), but would not invalidate the unbiasedness of the estimator. If, however, the choice of model form and the choice of predictors were sample-based, e.g., by choosing predictors by optimizing the predictive power of the model for the sample at hand, there would be a risk of overfitting and underreporting of uncertainty [72].
On this background, we found ourselves in a dilemma in this study. On one hand, we had no prior information about useful variables derived from the given RS data for AGB modelling for the particular forest types under study. Neither had we any experience with suitable model forms for the study area. On the other hand, if model selection and variable selection were optimized for the given sample, overfitting would be a likely consequence.
To balance these conflicting requirements, we first did a screening of the variables mentioned above to gain first-hand experience with the three types of satellite data for the current forest types. We then chose a model-form a priori, and allowed only a small number of predictors to be included in the model. In the modelling phase, we paid special attention to any sign of overfitting.
Thus, in the first phase of the analysis, Pearson's correlation coefficient was used to explore the relationships of individual independent variables with AGB. Those variables that had a significant correlation with AGB were used as potential variables for the AGB model fitting. Furthermore, correlation analysis was done for each pair of independent variables within each satellite data source to evaluate the level of intercorrelation between them. Results of the correlation analysis indicated that most of the variables were strongly intercorrelated ( Figure 2). Hence, variable screening was employed to reduce the redundant information emanating from those strongly intercorrelated variables. Results of the initial analysis using more complex models showed overfitting problems, which was manifested in precision difference between training and validation results for each model. Such severe overfitting was observed for models with more than two variables. Because of the risk of overfitting, we restricted the selection of independent variables in the models to a maximum of two variables only. The results from the analysis of models with more than two variables are not documented any further.
another one with the mean of SWIR1 SB were other candidate models of the S2 category.
From the PS data, the mean of G reflectance had the strongest correlation with AGB. Thus, one of the PS models contains independent variables composed of the mean reflectance of G SB and the standard deviation of the ASM texture variable of the NIR SB. The mean of B4ASM was the least intercorrelated with the mean of G SB. The other simple model was the model with a predictor variable of the mean of G SB reflectance only.  Table 4 for the descriptions of the notations used to represent the independent variables. Table 5 shows a detailed description of the candidate AGB models for each image type. Two candidate models were obtained from the L8 data. There was a marginal difference between the single and two-variable models with AIC of 1403.31 and 1402.68, respectively. The model calibration RMSE of the single and two-variable models were 70.22% and 71.06% of the mean AGB, respectively. Likewise, the respective model validation RMSE values were 73.23% and 73.31% of the mean AGB. As clearly revealed in these model metrics, there is concern of less responsiveness of the selected variables for the AGB estimates in the two-variable model. The presence of two variables in the model did not significantly improve the model performance. Therefore, the model with the mean of SWIR1 reflectance as the only predictor variable was selected for AGB estimation.

Selected AGB Models for Each Image Type
Three models were selected as candidates from the S2 variables. Two of them were with a single predictor variable while the third has two variables ( Table 5). The model with the mean of ExGI as a predictor variable had a larger validation RMSE (73.80%) than the other models. The model with the predictor variable of the mean of SWIR1 was better than the one with the mean of ExGI. However, the two-variable model had even greater performance among the S2 category of models. The twovariable S2 model with predictor variables of the mean of SWIR1 and standard deviation of GLI had the least AIC value among the models (1385.06) and minimal overfitting problem (Table 5). Additional indicators of the model fit and validation results of this model were better than the other models in the category. This model explained 40.96% of the variability in the ground reference AGB unlike the selected L8 and PS models, each of which explained less than 30%.
Two candidate models were obtained from the PS data. The two-variable PS model contains the mean of G reflectance and the standard deviation of B4ASM texture as predictor variables. However, this model revealed a severe overfitting problem. The model RMSE and validation RMSE were 70.19% and 79.48%, respectively. Thus, the single-variable model with the mean of G reflectance was selected for AGB estimation in this category. It had model calibration and validation RMSE of 71.79%  Table 4 for the descriptions of the notations used to represent the independent variables. Table 4. Correlation of relevant independent variables (see Tables 1-3 for definitions) derived from L8, S2 and PS images with AGB. The notations for SB and SI variables of all the image types are MMM_mean or MMM_std representing the mean and standard deviation of the variable MMM, in respective order. For texture variables of the PS data, BnXXX_mean is the mean of the mentioned (XXX) texture variable of the SB Bn where n = 1, 2, 3, 4 for B, G, R, NIR, respectively. Similarly, BnXXX_std is the standard deviation of the texture variable as described above for the BnXXX_mean, except replacing 'mean' by 'std'. VAR and MEA stand for the GLCM variance and mean, respectively. The relevant variables of each satellite data source were related to plot-level AGB using the logarithmic link function in a generalized linear model (GLM) of the form: where y i is ground reference AGB (Mg ha −1 ), β 0 is intercept, β i is the coefficient of the independent variable (X i ), and i is the index of an individual independent variable. This model form was chosen since it provides valid estimates where true zeroes are included in the estimate of AGB, which has positive continuous numerical values. A study of AGB prediction using topographic variables in human-impacted tropical dry forest landscapes of Mexico indicated that GLM estimation technique improved predictions [73]. Thus, the mean of SBs and SIs of L8 image were candidate independent variables for the L8 model. The mean and standard deviation of the SBs and SIs of the S2 image were candidate independent variables for the S2 model. The mean and standard deviation of SBs, SIs and texture features of PS bands were used as candidate independent variables for the PS model.

Model Validation
We evaluated the performance of the models using a leave-one-out-cross validation technique. The cross-validation was used to assess overfitting. Each model was validated in terms of coefficient of determination (R 2 ), root mean squared error (RMSE, %), mean deviation (MD, %), and Akaike Information Criterion (AIC) as determined by Equations (3)- (8). The AIC was used to evaluate the maximum likelihood of the model parameters. The maximum likelihood estimation enables choosing the parameter that makes the likelihood of having the observed data a maximum fit with the dependent variable (AGB) without causing an overfitting issue. When comparing models, the model with a smaller AIC is better than the one with a higher AIC.
where y i andŷ i are the ground reference and predicted AGB (Mg ha −1 ) in the ith SP; y is the mean of ground reference AGB ( Mg ha −1 ) of all SPs; n is the sample size; L β (k) is the likelihood function of the observations,β(k) is the maximum likelihood estimation of the parameter β given the number of parameters of k within the model. In addition to the validation metrics indicated above, we did qualitative evaluation based on a visual comparison between the predictions using the selected models in each satellite data source and false-color composite (i.e., band combination of NIR-R-G in the R-G-B channels) depiction of the S2 image.

Population-Level Estimation and Efficiency Assessment
Based on the SP inventory data, for the sample size of 111 plots of about 1000 m 2 area, the estimators of the mean AGB for the population and its variance were calculated by Equation (9) and Equation (10), respectively [71]:μ field = 1 n n i=1 y i (9) var(μ field ) = 1 n(n − 1) where y i is AGB (Mg ha −1 ) of the ith SP in the sample and n is the sample size. The 95% confidence interval (CI) ofμ field was calculated using Equation (11): where SE(μ field ) = v ar(μ field ) is the standard error (SE) ofμ field and t is student's t at a significance level of 0.05.
Similarly, we estimated the mean AGB for the entire study area using the selected regression model for each satellite data source. For this purpose, the study area was tessellated into grid cells of 31.64 × 31.64 m providing a total of N (141,604) population units. The size of the grid cells was chosen to be equivalent to that of the SPs. Area-weighted mean and standard deviation of the variables used in the regression models were extracted for each grid cell using QGIS. AGB was predicted for each population unit (i) in the map of the tessellated granules using the selected regression models for each satellite data source and is represented byŷ i . Because the prediction relied on field data collected based on probability sampling inside the population of interest, we adopted generalized model-assisted regression estimators. The mean and the variance estimates were computed using Equation (12) and Equation (13), respectively [71] (p. 231): where ε i and ε are the estimates of error at each data point (i) and the average, respectively. The SE of the mean AGB estimators (i.e., SE(μ field ) and SE μ image ) were calculated by taking the square root of the respective variance estimatorsv ar(μ field ) andv ar μ image .
The study assessed the gain in precision of AGB estimation with the use of the three types of RS data. The measure of quantifying such a gain in precision of using RS data over the pure field-based estimates was expressed using relative efficiency (REf). REf quantifies the magnitude of estimated variance of a remote sensing-assisted estimate of mean AGB to a field-based estimate. It was computed by Equation (14) as the ratio of the variance of the field-based estimates to the remote sensing-assisted estimates: When REf is greater than one, it is interpreted as the amount of additional precision gained due to the use of the RS data for estimating mean AGB.

Relationship of Independent Variables with AGB
Statistical test for significance of correlation coefficients of the relationship of individual RS variables with AGB demonstrated that many of the candidate variables were reasonably related to AGB. Correlation coefficients were translated to descriptors like 'weak', 'moderate' and 'strong' relationships according to the scheme used by [74]. The mean SB reflectance values of the three satellite data sources had negative moderate correlation with AGB (Table 4). On the other hand, the mean values of most SIs tend to show moderate positive relationships with AGB with some exceptions (for instance, mean ExGI of S2). It was revealed from the exploratory analysis that standard deviation of SIs of S2 and PS images and texture variables of PS images had moderate relationships with AGB.
For the L8 category of independent variables, the mean SIs were less correlated with AGB as compared to those of the SBs. Table 4 shows that from the SIs, the mean NDMI had considerable relationship while that of ARVI, NDVI and SR have weaker performance. The mean values of all SBs were moderately related to AGB with correlation coefficients ranging from −0.38 (NIR_mean) to −0.48 (SWIR1_mean). The strength of the association of AGB with NIR_mean was equivalent to that observed with the strongly correlated SI (i.e., the NDMI_mean), which was 0.39.
From the S2 variables, the mean of both SBs and SIs showed reasonable association with AGB. Similar to the L8 variables, the mean of SBs had a stronger relationship than that of the SIs except the mean ExGI, which had the strongest relationship. The peculiar behavior of ExGI comes from the fact that it is just a difference of SBs. Likewise; the standard deviation of the SIs (namely GLI, NDGI and VI) had strong positive associations with the dependent variable.
Similarly, for the PS variables, the mean of SBs of G, R and B showed the strongest relationship with AGB followed by that of VI and NDGI SIs. The mean and standard deviation of the texture data had moderate relationships with AGB (Table 4).

Variable Selection for the Prediction Models
Correlation analysis indicated that independent variables of each satellite data source were strongly intercorrelated ( Figure 2). Therefore, the variables that fit well with AGB in the GLM, and which had no significant collinearity problem, were selected for the AGB prediction models. As a result, the means of NDMI and NIR variables were less intercorrelated and became predictor variables for one of the L8 models. Besides, a simple model with the most strongly correlated variable (SWIR1) with AGB was considered as another candidate model in this category.
Similarly, the mean of SWIR1 and standard deviation of GLI were selected as predictor variables for the two-variable S2 model. The standard deviation of GLI had a strong positive correlation with AGB (Table 4) and was less correlated with the mean of SWIR1 variable (Figure 2), which was already in the model. Moreover, the single variable model with a predictor variable of mean of ExGI and another one with the mean of SWIR1 SB were other candidate models of the S2 category.
From the PS data, the mean of G reflectance had the strongest correlation with AGB. Thus, one of the PS models contains independent variables composed of the mean reflectance of G SB and the standard deviation of the ASM texture variable of the NIR SB. The mean of B4ASM was the least intercorrelated with the mean of G SB. The other simple model was the model with a predictor variable of the mean of G SB reflectance only. Table 5 shows a detailed description of the candidate AGB models for each image type. Two candidate models were obtained from the L8 data. There was a marginal difference between the single and two-variable models with AIC of 1403.31 and 1402.68, respectively. The model calibration RMSE of the single and two-variable models were 70.22% and 71.06% of the mean AGB, respectively. Likewise, the respective model validation RMSE values were 73.23% and 73.31% of the mean AGB. As clearly revealed in these model metrics, there is concern of less responsiveness of the selected variables for the AGB estimates in the two-variable model. The presence of two variables in the model did not significantly improve the model performance. Therefore, the model with the mean of SWIR1 reflectance as the only predictor variable was selected for AGB estimation.

Selected AGB Models for Each Image Type
Three models were selected as candidates from the S2 variables. Two of them were with a single predictor variable while the third has two variables ( Table 5). The model with the mean of ExGI as a predictor variable had a larger validation RMSE (73.80%) than the other models. The model with the predictor variable of the mean of SWIR1 was better than the one with the mean of ExGI. However, the two-variable model had even greater performance among the S2 category of models. The two-variable S2 model with predictor variables of the mean of SWIR1 and standard deviation of GLI had the least AIC value among the models (1385.06) and minimal overfitting problem (Table 5).
Additional indicators of the model fit and validation results of this model were better than the other models in the category. This model explained 40.96% of the variability in the ground reference AGB unlike the selected L8 and PS models, each of which explained less than 30%.
Two candidate models were obtained from the PS data. The two-variable PS model contains the mean of G reflectance and the standard deviation of B4ASM texture as predictor variables. However, this model revealed a severe overfitting problem. The model RMSE and validation RMSE were 70.19% and 79.48%, respectively. Thus, the single-variable model with the mean of G reflectance was selected for AGB estimation in this category. It had model calibration and validation RMSE of 71.79% and 75.17%, respectively. Overall performance of the selected PS model was slightly lower than the selected S2 model but similar to that of the selected L8 model ( Table 5).
The validation results of all the three selected models for AGB estimation indicated that the models have sensible performance in predicting AGB for the data with which they were not trained. The scatter plot of fitted versus ground reference AGB values shown in Figure 3 indicates a reasonable predictive power of the models given the complex settings of the study area. Pearson's correlation coefficients of the model predicted and ground reference AGB in the SPs revealed that the S2 model predictions were more correlated with the ground reference AGB than the other two models. The S2 model predictions had a correlation coefficient of 0.64 with the ground reference AGB.
The L8 and PS models had equivalent performance and explained a considerable amount of the variation in the FBSS estimate of mean AGB with R 2 of 29% and 27%, respectively, given the complex forest structure and topography in the study area.
Although the general trend of the error distribution of the three selected models looks similar, prediction errors of the L8 and PS models spread out at small and large AGB more than the S2 model did (Figure 3). The extents of deviation of predicted values from the ground reference AGB differ for each model especially at the smaller and larger AGB values. With this variability maintained, the selected models of all the three-image data inflated predictions of small AGB, particularly those below approximately 300 Mg ha −1 . For SPs with large AGB, the predictions using all the three models were smaller than the ground reference values.
The L8 and S2 models had smaller prediction error at the small AGB end than at the large AGB levels. Generally, the predictive power of the S2 model prevailed over that of the other models. Table 6 shows the estimated mean AGB, estimates of mean deviation, SE of the mean AGB estimates and REf for the selected models of the three image categories presented in Table 5. The estimates of mean AGB were 179.67 Mg ha −1 , 177.79 Mg ha −1 and 184.27 Mg ha −1 when using the L8, S2 and PS model predictions, respectively. The model-assisted estimates of the mean AGB for all the three categories of models were within 95% CI of the mean AGB estimate based on the field data only (i.e., 155.15-213.76 Mg ha −1 ). The estimated mean AGB using the PS model was closer to the field-estimated mean AGB (i.e., 184.35 Mg ha −1 ) than the estimates using the other models. The estimated mean AGB using the L8 and PS models had the largest and smallest MDs, respectively. The AGB estimate based on the PS model was relatively less precise followed by the L8 model. The estimation results revealed that the L8 and PS models resulted in equivalent estimation efficiencies (i.e., 1.40 and 1.37, respectively). and 75.17%, respectively. Overall performance of the selected PS model was slightly lower than the selected S2 model but similar to that of the selected L8 model ( Table 5).

Estimation and Mapping of AGB Using the Selected Models
The validation results of all the three selected models for AGB estimation indicated that the models have sensible performance in predicting AGB for the data with which they were not trained. The scatter plot of fitted versus ground reference AGB values shown in Figure 3 indicates a reasonable predictive power of the models given the complex settings of the study area. Pearson's correlation coefficients of the model predicted and ground reference AGB in the SPs revealed that the S2 model predictions were more correlated with the ground reference AGB than the other two models. The S2 model predictions had a correlation coefficient of 0.64 with the ground reference AGB.
The L8 and PS models had equivalent performance and explained a considerable amount of the variation in the FBSS estimate of mean AGB with R 2 of 29% and 27%, respectively, given the complex forest structure and topography in the study area. Although the general trend of the error distribution of the three selected models looks similar, prediction errors of the L8 and PS models spread out at small and large AGB more than the S2 model did (Figure 3). The extents of deviation of predicted values from the ground reference AGB differ for each model especially at the smaller and larger AGB values. With this variability maintained, the selected models of all the three-image data inflated predictions of small AGB, particularly those below approximately 300 Mg ha -1 . For SPs with large AGB, the predictions using all the three models were smaller than the ground reference values.
The L8 and S2 models had smaller prediction error at the small AGB end than at the large AGB levels. Generally, the predictive power of the S2 model prevailed over that of the other models.  Table 6. Estimated mean AGB (Mg ha −1 ), mean deviation (MD) in Mg ha −1 , standard error (SE) of the mean estimates (Mg ha −1 ) and relative efficiency (REf) when using the selected models to assist in the estimation. The estimate based on the S2 model was the most precise among the three model-assisted AGB estimates with SE of 11.40 Mg ha −1 . As a result, the REf of the mean AGB estimate using the S2 model (i.e., 1.68) was greater than what we obtained by using the other two models.

Estimator Data Source Estimated Mean AGB Estimated MD SE Ref
Visual inspection of the predicted AGB using the three selected models and the false-color composite of the S2 image shows convincing AGB distribution across the landscape. As expected, the patches of bare land (non-forest areas) in the study area (shown in different shades of grey in Figure 4D) have small AGB predictions using all the models (shades of yellow in Figure 4A-C) while the dense forest areas (colored red in Figure 4D) yielded greater predicted AGB values (shades of green in Figure 4A-C). The map revealed that AGB predictions using the selected models of the three satellite image types had many similarities, which also was confirmed by similarities in the estimated uncertainties ( Table 6). The distribution patterns of AGB predictions in the maps ( Figure 4A-C) indicated the spatial consistency of AGB predictions across the area for the selected models.

Variable Exploration for Estimating AGB and Model Selection
The observed moderate relationships of independent variables of the RS data with AGB demonstrated the potential of optical RS data for developing models to enhance AGB estimation. The

Variable Exploration for Estimating AGB and Model Selection
The observed moderate relationships of independent variables of the RS data with AGB demonstrated the potential of optical RS data for developing models to enhance AGB estimation. The observed negative correlation coefficients between the mean of SB reflectance and AGB agree with results of similar studies conducted in various forest types [19,75,76]. The negative correlation coefficients indicate the inverse relationship between reflectance values from the SBs and AGB. This relationship in the current study could be explained by a shadow effect within the complex forest stands where AGB is large [77,78]. The presence of scattered big trees in SPs with large AGB results in large shadows. Additionally, such an effect might be related to large canopy water content, which is directly linked to photosynthetic efficiency [79]. The reflectance of the SBs from uniform forest stands like young plantations is large but they have relatively small AGB.
The positive relationship of most of the SIs with AGB found in this study is in accordance with previous research findings [19,34,76]. Besides the mean of SIs, the standard deviation of some SIs had also remarkable potential to relate with AGB. Naesset et al. [16] got similar results in Tanzania. Reviewed literature indicated that application of some of the SIs like the GLI, ExGI and NDGI have been limited to assessing grass biomass and crop cover or yield estimation. However, the current study showed that they had great potential to predict AGB in this type of forest. Thus, an in-depth study is required to understand the potential of such SIs for AGB estimation in different forest types.
It was revealed from the correlation analysis that most of the predictor variables in each satellite data type were intercorrelated. Lu et al. [80] found a similar result for estimation of AGB in wheat using an unmanned aerial vehicle. Among the different SBs the visible and SWIR bands, which are affected by atmospheric interference and shadow, were more strongly intercorrelated [78]. Besides, the SIs and texture variables were derived from these interrelated SBs. Therefore, the observed intercorrelation among the independent variables was likely to happen. This suggests the importance of a careful screening of RS variables for AGB modelling.
Furthermore, inter-resolution comparison of SBs showed that the limited spectral properties of the PS images might have restricted their potential to characterize AGB. For example, AGB correlated similarly with the G SB from each of the three image sources regardless of the differences in spatial resolution. The study results showed that the same SB across the resolution gradient characterized AGB similarly, indicating only a minor impact of pixel resolution on the quality of the AGB models if only SBs are related to AGB (Table 4). We observed that the S2 data contain a range of SBs that were more sensitive to AGB than the PS data, which have a higher spatial resolution.
Based on the relationship of the independent variables with AGB, we identified useful variables and models for each satellite data source. For example, exploration of the L8 data revealed that the NDMI showed a stronger correlation with AGB than other SI variables including the NDVI. This might be due to the improvements in the NDMI to detect leaf water content at the canopy level [7]. Previous research indicated that NDMI is useful for predicting forest attributes, including biomass [35,81]. The NDVI, which is the most popular SI for AGB modelling mostly in the temperate and boreal forests, was not a good predictor of AGB in the current study. Sader et al. [32] got similar results indicating unsuitability of NDVI for estimating AGB in tropical dense forests.
However, for the L8 data, AGB had a stronger correlation with SBs than the SIs mentioned above. Even among the SBs, AGB strongly correlated with the mean of SWIR1. The significance of the SWIR1 variable for AGB modelling was according to the results of other studies [23]. A study of biomass estimation using RS data in India indicated that biomass models using the SWIR bands were more reliable than those using short-wavelength SBs like the visible bands, which are more sensitive to atmospheric effects [82,83]. For green vegetation, reflectance in the SWIR spectral regions is controlled by the amount of water in the leaf biomass of the canopy. There is low diffuse of light at the SWIR wavelengths, and hence shadows are contrasted. The presence of thick layers of fragmented tree canopy and shadows in SPs with large AGB yielded low reflectance in the SBs, including the SWIR1 SB as indicated with the negative coefficient in the selected model.
The S2 model used for AGB estimation was the two-variable model with the independent variables of the mean of SWIR1 and standard deviation of GLI. Inclusion of the standard deviation of the GLI variable in the S2 model indicates the ability to capture spatial variation in canopy structure in the forest as the GLI can identify green leaves and stems from the background soil surface [60]. This variable may reflect the level of disturbance, terrain variation or presence of very big scattered trees in the natural forest. This variable signifies the importance of using measures of variability derived from higher resolution images in AGB modelling. The mean of ExGI was also another variable from the S2 data sensitive to AGB variability.
Besides the SBs and SIs, texture variables of the PS images had considerable potential for AGB modelling. The standard deviation of B4ASM that was included in one of the PS models, indicating the importance of high-resolution images for AGB modelling. Image texture variables like the ASM, describe the spatial arrangement of pixels with varying intensities that resulted in different AGB. The texture variables were able to differentiate between heterogeneous and homogenous surfaces, which prevailed in the disturbed natural forest patches and young plantation forests, respectively. This might be the reason for the observed strong positive relationship of the standard deviation of B4ASM with AGB. Improvement in the model performance by including this texture variable was in line with the findings of other studies [24,39,84]. Nevertheless, the two-variable model containing the standard deviation of B4ASM was subject to overfitting as compared to the reduced model with the mean of the G reflectance variable (Table 5).

Model Characteristics and Their Contribution to Enhance AGB Estimation
Generally, the selected L8 and PS models explained some proportion of the variability in the ground reference AGB that was better than the results from [16] although the nature of the current forest and terrain configuration was complex. The calibration RMSEs of the L8 and PS models were 71.06% and 71.79% of the mean of the ground reference AGB, respectively. This was comparable with the results of other studies conducted even in intermediate vegetation cover conditions where it is easier to get a stronger relationship between image data and AGB [70]. The REf of the AGB estimates based on the selected L8 and PS models were 1.40 and 1.37, respectively. That means the selected L8 and PS models could reduce sample sizes to 71% and 73%, respectively, of the field sample size to get the same precision with the FBSS estimates. These amounts (or proportions) of the variability in the field estimates remained unexplained when the model-assisted estimation was applied. Consequently, there was a similarity in improving the mean AGB estimates based on the L8 and PS models. The REf when using the L8 model in this study was slightly larger than the findings by Naesset et al. [16] for Miombo woodlands in Tanzania using the global Landsat products. The selected S2 model contributed more strongly to improve the precision of the AGB estimates than the L8 and PS models. This improvement in estimation efficiency contributes to reducing the number of field SPs required to attain the same precision, to approximately 59% of the sample size required for a pure FBSS estimate. The REf when using the S2 model was smaller than that of the RapidEye images used for AGB estimation in the Miombo woodlands in Tanzania [16]. This might be attributed to the heterogeneity of the forest in the current study or the interaction effect of forest types and spatial resolution of the images. Besides, [16] stated that the small study area covered in their study might have resulted in overly optimistic results because the RS data were very homogenous since they came from only a single scene. On the other hand, the results of the current study were similar to the findings by Navarro et al. [21] who studied AGB of mangrove plantations using S2 images in Senegal. Thus, the findings of the current study are reasonable given the heterogeneity of the terrain and forest conditions, which influence the relationship of image data and AGB [83].
Although there might be some variations between the models in this regard, they were able to predict only to a limited range of the ground reference AGB. This shows a saturation problem for which canopy shadow is mainly responsible in the SPs with large AGB. Similar studies in primary and successional forests in Brazil indicated that shadows were among the main factors resulting in data saturation, particularly in natural forests with large AGB [77]. Modelling of tree heights might also contribute to the low/moderate estimation efficiencies in these models. There is uncertainty inherent in the field measurements of AGB, which could be caused by the AGB estimation procedure (e.g., errors in measurement and in the allometric models). The error in the height-DBH models could increase the uncertainties and limit performance of the models. Therefore, future efforts should focus on synchronizing other auxiliary variables like canopy density and canopy height from airborne laser scanning data with the identified variables to improve the model performances.
During the fieldwork, understory vegetation was observed in SPs (see Figure 5). The inflated predictions at SPs with small AGB might be attributed to this phenomenon. The field inventory was limited to trees with DBH ≥ 5 cm and did not account for the understory vegetation. Dense understory vegetation composed of saplings, shrubs, lianas and herbaceous species covered most of these SPs and challenged our movement during the fieldwork. The biomass in the understory vegetation, which was not accounted for in the ground reference AGB, could have had a major influence on the SB reflectance values and hence in all the RS variables. This might partly explain the moderate improvement in the precision of the model-assisted estimates of AGB compared to the pure field-based estimate. As shown in Table 6, the PS model-assisted estimates had negative MD indicating the greater effect of the inflated predictions at SPs with small AGB than the reduced predictions at the SPs with large AGB. Thus, the effect of understory vegetation on the relationship between image-derived variables and AGB was more obvious when using the high-resolution PS images. A greater compliance of the RS data with AGB would happen for homogeneous forests in which the understory vegetation cover is minimal and the forest canopy cover is uniform. Therefore, further studies are needed in pure plantation forests to attain an optimum efficiency of RS data for AGB estimation beyond the ones we got in this study.  Generally, the findings of the current study were encouraging. We identified relevant variables extracted from RS data for AGB estimation. The selected models of each satellite data source based on the identified variables provided reasonable improvements in AGB estimations, which were reinforced by other research findings. The freely available S2 data were particularly useful. The research results revealed that S2 images possess sensible spectral and spatial properties for AGB estimation. The results of this study will help to satisfy the existing demand for forest carbon stock assessment by the national REDD+ program in Ethiopia. Enhanced forest information using the freely available data sources like the S2 would help to improve sustainable forest management and encourage results-based payments for those who properly manage their forest resources according to established principles like the REDD+ schemes. Generally, the findings of the current study were encouraging. We identified relevant variables extracted from RS data for AGB estimation. The selected models of each satellite data source based on the identified variables provided reasonable improvements in AGB estimations, which were reinforced by other research findings. The freely available S2 data were particularly useful. The research results revealed that S2 images possess sensible spectral and spatial properties for AGB estimation. The results of this study will help to satisfy the existing demand for forest carbon stock assessment by the national REDD+ program in Ethiopia. Enhanced forest information using the freely available data sources like the S2 would help to improve sustainable forest management and encourage results-based payments for those who properly manage their forest resources according to established principles like the REDD+ schemes.

Conclusions
Optical RS images from L8, S2 and PS satellites were studied to identify relevant RS predictor variables that could be used to enhance AGB estimation in a dry Afromontane forest. Most of the SBs, some SIs and texture variables (listed in Table 4) were found to be promising variables for predicting AGB. Although some of them were not selected in the models used for assisting AGB estimation, we identified variables including the mean of GLI, ExGI and NDGI that were seldom used for AGB modelling but are highly correlated with AGB. We recommend a detailed investigation of the importance of these variables for AGB assessment in various forest conditions.
The simple models selected for each satellite data source enhanced AGB estimation. Of the variables used in the models, the SWIR1 SB, which lacks in the PS data, was a useful variable of the L8 and S2 images for AGB estimation in this forest type despite the huge differences in pixel resolution among the image types. The study suggested that the additional spectral information of L8 and S2 images was more determinant of AGB estimation than the small pixel size of the PS images.
The use of RS data for AGB estimation improved the precision of estimates. Thus, the remote sensing-assisted estimation techniques used in this study will complement the FBSS estimates of AGB by improving precision. The model-assisted estimation will reduce sample sizes to obtain a similar estimation efficiency with the field survey. However, the models used for AGB estimation in this study revealed saturation problem. Therefore, future studies should focus on refining these limitations using a synergy of different data sources to enhance the estimation efficiency of AGB models beyond the ones achieved in the current study.
The methods used in this study could be adopted to similar conditions in forests that have limited application of RS data. The potential predictor variables derived from optical satellite images for biomass estimation were identified from studies showing global experiences. Exploratory data analysis was used to identify relevant predictor variables for biomass estimation in the current study site. Choice of a model form that is important for biomass required understanding the characteristics of data types. The selected models for each image type predicted biomass with estimation efficiencies comparable with those obtained in other forest types. These methods contain a unique mix of techniques capable of using satellite images for biomass estimation in a data scarce forest type.