A Novel Method for Estimating Spatial Distribution of Forest Above-Ground Biomass Based on Multispectral Fusion Data and Ensemble Learning Algorithm

: Optical remote sensing technology has been widely used in forest resources inventory. Due to the inﬂuence of satellite orbits, sensor parameters, sensor errors, and atmospheric effects, there are great differences in vegetation spectral information captured by different satellite sensor images. Spectral fusion technology can couple the advantages of different multispectral sensor images to produce new multispectral data with high spatial and spectral resolution, it has great potential for improving the spectral sensitivity of forest vegetation and alleviating the spectral saturation. However, how to quickly and effectively select the multi-spectral fusion data suitable for forest above-ground biomass (AGB) estimation is a very critical issue. This study proposes a scheme (RF-S) to comprehensively evaluate multispectral fused images and develop the appropriate model for forest AGB estimation, on the basis of random forest (RF) and the stacking ensemble algorithm. First, four classic fusion methods are used to fuse the preprocessed GaoFen-2 (GF-2) multispectral image with Sentinel-2 image to generate 12 fused Sentinel-like images. Secondly, we apply a comprehensive evaluation method to quickly select the optimal fused image for the follow-up research. Subsequently, two feature combination optimization methods are used to select feature variables from the three feature sets. Finally, the stacking ensemble algorithm based on model dynamic integration and hyperparameter automatic optimization, as well as some classic machine learners, are used to construct the forest AGB estimation model. The results show that the fused image NND_B3 (based on nearest neighbor diffusion pan sharpening method and Band3_Red) selected by the evaluation method proposed in this study has the best performance in AGB estimation. Using the stacking ensemble method and NND_B3 image, we get the highest estimation accuracy, with the adjusted R 2 and relative root mean square error ( RMSEr ) of 0.6306 and 15.53%, respectively. The AGB estimation RMSEr of NND_B3 is 19.95% and 24.90% lower than those of GF-2 and Sentinel-2, respectively. We also found that the multi-window texture factor has better performance in the area with low AGB, and it can suppress the overestimation signiﬁcantly. The AGB spatial distribution estimated using the NND_B3 image matches the ﬁeld observations well, indicating that the multispectral fusion image combined with the Stacking algorithm can increase the accuracy and saturation of the AGB estimates.


Introduction
The forest ecosystem is very important to the Earth's ecosystem. Its biomass and carbon storage play a very important role in global climate change and material circulation, the image spectral information [24]. However, only the Gram Schmidt fusion algorithm is used in their study, and there is no comparative analysis about other fusion methods or other optical images. Therefore, it is necessary to explore the utility of other classic fusion algorithms for fusing Sentinel-2 and GF-2 multispectral images for forest AGB estimation. Moreover, effective evaluation and accurate selection of multi-spectral fusion data suitable for forest AGB estimation is a key issue.
Feature selection is the key to improve the performance of the forest AGB estimation model [46][47][48][49][50][51][52]. However, there is little evidence to indicate that the features screened by a certain method have a good accuracy of AGB or GSV estimation for all regression models [24,52]. Li et al. [24] proposed a feature variable screening and combination optimization procedure based on the distance correlation coefficient and k-nearest neighbor algorithm (DC-FSCK). This algorithm considers the correlation, heterogeneity and combination optimization characteristics between feature variables, and can select the best feature combination for the forest GSV estimation research. The DC-FSCK method has achieved very satisfactory results when the k-nearest neighbor (KNN) algorithm is used [53]. However, when the random forest (RF) regression algorithm is used, the performance of DC-FSCK is not ideal. Therefore, in order to improve the robustness of the feature variable combination method, we need to optimize the regression model of the algorithm based on the original DC-FSCK algorithm framework in order for it to adapt more application scenarios.
The commonly used models for AGB estimation include parametric regression models and non-parametric machine learning regression models [11,20]. Parameter regression models, such as linear regression and perceptron, estimate AGB by constructing a regression formula on the basis of the relationship between the measured AGB and remote sensing feature variables, topographic factors, and forest stand parameters. Non-parametric machine learning models, such as RF regression and support vector regression (SVR), can obtain different fitting functions by the training data, which can finally predict the target value [23,[43][44][45][46][47][48][49][50]. Ensemble machine learning algorithms learn training samples by constructing and integrating multiple learners, which have higher accuracy than traditional non-parametric or parametric methods in the forest parameter estimation based on small sample data [24,43,[51][52][53][54]. Generally, there are three ensemble methods: stacking, bagging (e.g., RF), and boosting (such as adaptive boosting, gradient boosting decision tree, extreme gradient boosting, and categorical boosting) [52]. The stacking algorithm has good performance in forest GSV prediction and vegetation classification research [24,[55][56][57], but there is no report on its application in forest AGB estimation.
In short, remote sensing data sources, feature variable combinations, and the estimation models influence estimating the forest AGB using remote sensing data [20,36,44]. Therefore, this study will solve these problems by applying the following steps.
(1) GF-2 and Sentinel-2 multispectral images are fused by four classic fusion algorithms to get the Sentinel-like images. Then, these fused images are assessed by a comprehensive evaluation method, which using the image information entropy, grayscale mean, standard deviation, average gradient, and image-based model cross-validation estimation error as the comprehensive evaluation index. Hence, the fused image suitable for AGB estimation is screened out.
(2) Feature variables are extracted from the fused image and terrain feature derivation factors (e.g., elevation, slope, aspect), and two feature combination optimization methods are used to screen feature variables for AGB estimation.
(3) The ensemble machine learning algorithm is utilized to build a forest AGB estimation model based on the selected feature combination, and is compared the performance with other machine learners.
We expect that the combination of fused images, new feature selection method and ensemble machine learning algorithm will yield a quickly and highly accurate AGB map of forests.

Study Area
This study was conducted in a state-owned forest farm (Huangfengqiao) located in the southeast of Hunan province, China ( Figure 1). It covers the middle of Luoxiao Mountains and the southwest of Wugong Mountain. There are many low mountains with the elevation varying between 1270 m and 115 m (Figure 2a). It has a subtropical monsoon humid climate, with an average annual temperature of 17.8 • C, an annual precipitation of 1410.8 mm, and an annual frost-free period of about 292 days. The dominant tree species is Chinese fir (Cunninghamia Lanceolata), and there are also Betulaceae, Camphor (Cinnamomum camphora), etc. (Figure 2b) [29]. The forest farm has a forest GSV about 891,000 cubic meters, and a forest coverage rate of 90.7% [29].

Field Plot Data Collection
Chinese fir plantation is distributed in the north, east and south of the study area. A total of 50 plots of Chinese fir were measured by field investigation using a random stratification sampling from 2016 to 2017. Due to the scarcity and inaccessibility of woodlands above 800 m, all the selected plots are below 800 m ( Figure 2a). The plots are 20 m × 20 m or 30 m × 30 m large, depending on the topographic features and tree stand density ( Figure A1a). The Zenith15A real-time kinematic (RTK) system was used to receive the signals from Hunan Satellite Navigation and Positioning Public Service Platform (HNCORS) and work together with the total station ZT20 to accurately measure the positions of the four corner points of the sample plot to determine the plot boundary ( Figure A1b,c)). In each plot, the information of all standing living trees with the diameter at breast height (DBH) no smaller than 5 cm were measured and recorded, including tree DBH, height, and topographic factors (e.g., slope, aspect). As shown in Table 1, the tree AGB value can be obtained by DBH and tree height [58]. The values of the AGB measurements at all plots are shown in Table 2.

Satellite Image Collection and Pre-Processing
We collected the L2A level product (atmospheric correction) of one Sentinel-2 image obtained on February 14, 2017 (https://scihub.copernicus.eu/, accessed on 17 March 2019), and six GF-2 images dated on 8 December, 2016 for this study (http://www.cresda.com/ CN/, accessed on 10 May 2019). The GF-2 satellite is the first optical remote sensing satellite with a spatial resolution of finer than 1 m independently developed by China. It is equipped with two high-resolution cameras of 1-m panchromatic and 4-m multispectral (blue, green, red, and near infrared) images. It has greatly improved the satellite's comprehensive observation efficiency, which has reached the international advanced level [59]. The Sentinel-2 satellite carries a multispectral imager (MSI), which can cover 13 spectral bands with the ground resolution of 10 m, 20 m, and 60 m, including three vegetation red edge bands and three short wave infrared bands that may improve forest AGB estimation accuracy [59]. The DEM data with the spatial resolution of 30 m × 30 m are obtained (http://www.gscloud.cn/, accessed on 15 June 2020) for terrain correction of optical satellite images.
Remote Sens. 2021, 13, 3910 5 of 28 tains and the southwest of Wugong Mountain. There are many low mountains with the elevation varying between 1270 m and 115 m (Figure 2a). It has a subtropical monsoon humid climate, with an average annual temperature of 17.8 °C, an annual precipitation of 1410.8 mm, and an annual frost-free period of about 292 days. The dominant tree species is Chinese fir (Cunninghamia Lanceolata), and there are also Betulaceae, Camphor (Cinnamomum camphora), etc. (Figure 2b) [29]. The forest farm has a forest GSV about 891,000 cubic meters, and a forest coverage rate of 90.7% [29].

Data Preparation
Chinese fir plantation is distributed in the north, east and south of the study area. A total of 50 plots of Chinese fir were measured by field investigation using a random stratification sampling from 2016 to 2017. Due to the scarcity and inaccessibility of woodlands above 800 m, all the selected plots are below 800 m (Figure 2a). The plots are 20 m × 20 m

The RF-S Model
In order to improve the accuracy and to alleviate saturation problem of forest AGB estimation, this research develops a novel integrated scheme on the basis of RF and stacking integration algorithm (referred to as the RF-S model hereafter). As shown in Figure 3, the RF-S model includes four steps: (1) Gram Schmidt (GS), Nearest Neighbor Diffusion pan sharpening (NND), Wavelet Resolution Merge (WRM), and Brovey Transform (BT) are applied to fuse each spectral band (Bule, Green and Red) of GF-2 with the Sentinel-2 image; (2) Predictor variables (Table 3) of each image are selected by the RF method according to the importance, and the random forest regression (RFR) algorithm is used to build the AGB estimation model and obtain the relative root mean square error (RMSEr), then a comprehensive evaluation index is employed to assess all fused images to quickly select the optimal image for further processing; (3) The proposed feature variables combinatorial optimization method that is based on KNN and RFR algorithm is used to choose the best feature variables from the optimal image; (4) The stacking algorithm is utilized to build the AGB estimation model and map the AGB distribution of the study area.
In this study, we compare four image fusion algorithms, three feature variable sets, two feature selection methods, and four AGB estimation models, and get the best solution for AGB estimation of Chinese fir plantation.

Multispectral Image Data Fusion
In this study, in order to improve the resolution of the multi-spectral image, we used the GS method to fuse the GF2 panchromatic image with the multi-spectral (blue, green, red) images. Then, in order to couple the two sensors data and increase the image information, the fused GF-2 multispectral images were fused with the Sentinel-2 images to generate Sentinel-like multispectral images with a spatial resolution of 1 m. Compared with the original GF-2 and Sentinel-2 images, the fused Sentinel-like images contain more details and more spectral information. In order to find a multi-spectral image fusion method suitable for forest AGB estimation, we compare four image fusion algorithms, GS, NND, WRM, and BT. In this study, each multispectral band (B1_blue, B2_green, and B3_red) of the GF-2 image was fused with the Sentinel-2 image, and the obtained images contain 10 bands including 4 vegetation red edges ( Table 4). The fused images are denoted by the fusion method and band name. For example, using the GS method to fuse the B1_blue image of GF-2 with Sentinel-2, the obtained image is denoted as GS_B1. We assume that the spectral information contained in each pixel in the image has a specific correlation with the average AGB of trees per unit area, and this information will not be lost after the image is resampled. Therefore, after multispectral data fusion, this study needs to resample the fused high-resolution images to a resolution similar to the sample plot size. Then, the regression model is established by using the measured AGB of sample plots and pixel spectral information.

Vegetation Indices Equation Reference
Normalized difference vegetation index NDV I = Band N IR −Band RED Band N IR +Band RED [58] Similar normalized difference vegetation indices Simple two-band ratios RV I i_j = Band i Band j [58] Enhanced vegetation index Difference vegetation indices DV I i_j = Band i − Band j [23] Soil adjusted vegetation indices Band N IR +Band RED +k [24] Atmospherically resistant vegetation index Modified simple ratio

The Fused Image Feature Extraction
In this study, the spectral features of different remote sensing images are extracted, and the vegetation indices are calculated ( Table 3). The method of "extracting multiple values to points" in the ArcGIS 10.2.2 software is used to interpolate and extract the feature variables of the spectrum and vegetation index of each forest plot. Finally, the extracted variables are combined with the measured AGB data to generate a sample dataset.

The Feature Selection and the AGB Estimation RMSEr Calculation for Each Fused Image
Using the sample data set and RF regression algorithm, we establish the AGB estimation model that contains k decision trees. We use the out-of-bag (OOB) data to calculate the OOB data error of each decision tree in the RF regression model (Equation (1)), denoted as ErrOOB1 i . Then, we traverse all the feature variables, and randomly add noise to the feature variable V i , and calculate the OOB data error again, denoted as ErrOOB2 i . In this way, k decision trees can get k ErrOOB1 and ErrOOB2. As shown in Equation (2), the importance of the feature variable V i can be described as I M_RFerr V i by calculating the magnitude of the error change before and after adding noise [34]. This is called the RF mean decrease accuracy method and denoted as RF_MDA. The advantage of the RF_MDA method is that it can quickly and accurately measure the importance of feature variables. Obviously, for unimportant variables, disrupting the original order of variables will not Remote Sens. 2021, 13, 3910 9 of 28 have much impact on the accuracy of model estimation, but for important variables, this will significantly reduce the accuracy of model estimation.
where n, is the number of samples, y i andŷ i , are the observed and estimated AGB, respectively.
where p, is the number of feature variables. In order to quickly select the most suitable feature variables for forest AGB estimation, all feature variables of each image are ranked by the importance index of the RF method. The most important features are selected in turn to form various feature variable combinations. For each feature combination, we built a model by the RFR algorithm to predict the AGB. The RMSEr between the predicted and observed AGB is calculated by the leave one out cross validation (LOOCV). Each iteration of LOOCV method leaves only one sample as the test set and other samples as the training set. If there are k samples, the method needs to train K times and test K times, so as to maximize the use of all samples.
where y, is the mean of observed AGB values of all sample plots.

Image Evaluation and Selection
The selection of remote sensing images is very important for the estimation of forest AGB. This research used the image information entropy, grayscale mean (Mean), standard deviation (SD), average gradient (AG), and image-based model cross-validation RMSEr as the comprehensive evaluation index to select the optimal images for forest AGB estimation.
The information amount increase is an important factor in evaluating the fusion effect, which can be calculated by information entropy as follows: We also analyze whether the fused image has more spatial details and texture information than the original image. Generally, image brightness can be quantified by indicator of grayscale mean, the greater the mean value, the better the image brightness. The standard deviation can be used to evaluate the gray dispersion of the image. The greater the standard deviation value, the greater the image contrast. The average gradient of the image can reflect the definition of the image to a certain extent. The larger the average gradient of the image, the more spatial details will be reflected [42].

Forest AGB Estimation Modeling Based on the Selected Optimal Fused Images
Texture feature variables are very helpful for remote sensing data modeling [20,24], so we extract texture features from the selected fused images and terrain factors derived from DEM data, including elevation, slope, aspect, and topographic wetness index (TWI) [62]. Relevant studies show that macro topographical factors (e.g., TWI) are related to regional forest AGB [63]. TWI is a physical index of the impact of regional topography on runoff flow direction and accumulation, which is helpful to identify rainfall runoff patterns, potential areas of increased soil water content, and ponding areas. Generally, when other conditions of forest (e.g., environmental and climatic) are the same, the larger TWI is more conducive to the growth of trees.
where α, is catchment area per unit contour length, β is the steepest outward slope of each pixel. These features combine with the vegetation indices and the measured forest plot AGB value to form a training sample data set. Then, we select the optimal combination of feature variables for AGB estimation.

Feature Variable Extraction
In this study, the gray level co-occurrence matrix (GLCM) is used to measure texture of the optical images and terrain data [20]. By calculating the mean, variance, homogeneity, contrast, dissimilarity, entropy, second moment, correlation, we obtain the direction, field and change range of image gray, which reflect the correlation between texture gray levels well. We extract the textural images using the GLCM with step size [1,1] and window size (3 × 3, 5 × 5, 7 × 7, 9 × 9). In order to analyze the contributions of vegetation indices, single-window texture factor and multi-window texture factor to the forest parameter estimations, we design three sets of feature variables, F1, F2, and F3 (Table 4).

Feature Variable Combinations
The DC-FSCK method was proposed by Li et al. to choose the optimal feature variable combination for GSV estimation [24], but the feature variables selected by one method may not bring good accuracy of AGB estimation in all regression models [24,52]. In addition, RFR algorithm has excellent performance in many applications of forest mapping [34]. Li et al. [58] compared the performance of six feature variable selection methods in forest GSV estimation, and found that the feature variable combination optimization method based on RFR model performs best. In order to improve the robustness of the feature variable selection methods, we replace the KNN in DC-FSCK method with the RFR algorithms ( Figure A2). Therefore, in the second stage of the RF-S model, for each data scenario, we apply the KNN-based method and the RFR-based method to screen feature variables for forest AGB estimation.

Stacking Ensemble Algorithm
The estimation result of a single model has the disadvantages of one-sidedness and contingency. Therefore, the ensemble machine learning model that integrates the results of multiple models for prediction has better performance than a single model. Generally, there are three types of ensemble machine learning regression algorithms, including bagging, boosting, and stacking [59]. Stacking generalization algorithm was first proposed by Wolpert in 1992, and he believes that it is similar to cross validation, which is integrated through "winner takes all" [57]. Bagging and boosting algorithms usually took the decision regression tree as the basic model [56,57]. So they are the integration of similar models, while stacking algorithms are relatively more flexible, which can be the ensemble of similar models or heterogeneous models. This research plans to adopt the integration algorithm of the stacking machine learning model that combines the prediction results of the basic model to realize accurate AGB estimation. The base model, meta-model, and hyperparameter optimization are the keys for the stacking integration algorithm. In this study, the base models are SVR, RFR, and KNN, the meta-model is least absolute shrinkage and selection operator (LASSO). Therefore, this study comprehensively considers the coupling and complementarity between these models, and realizes the integration optimization of the regression model through the iterative selection of models and the automatic adjustment of hyperparameters (Figure 3 Stage2).

Model Evaluation and Application
In this study, four methods including KNN, SVR, RF, and the stacking algorithm are employed to conduct AGB modeling and estimation using NND_B3, GF-2, and Sentinel-2 image data. These results are assessed by the LOOCV method using the indexes of coefficient of determination (R 2 ), adjusted R 2 , Pearson correlation coefficient (r), RMSE, RMSEr, and mean absolute error (MAE).
where n and p, are the number of samples and feature variables, respectively.ŷ i , is the mean of estimated AGB values of all sample plots. Finally, the data scenario and estimation model with larger adjusted R 2 and smaller MAE and RMSEr are selected to map the Chinese fir plantation AGB in the study area.

Twelve Sentinel-Like Images Generated by Four Fusion Methods
In this study, four classic pixel-level fusion algorithms (GS, NND, WRM, and BT) are used to fuse three GF-2 multispectral images with Sentinel-2 images. Compared with the original Sentinel-2 image ( Figure A3), all Sentinel-like RGB true color images ( Figure 4) have enhanced clarity. Specifically, the GS and NND fusion images (Figure 4(a1-a3,b1-b3)) have more obvious texture details, better clarity. The overall tone of the WRM fusion images (Figure 4(c1-c3)) is natural, close to Sentinel-2 image, indicating high spectral retention. However, these images have no obvious texture details and have a lot of noise and outliers. The BT fusion images (Figure 4(d1-d3)) have more obvious texture details, but lower tone and larger spectral distortion than the original Sentinel-2 image.
As shown in Figure 5, Band 8 (Vegetation Red Edge 4) of all the Sentinel-like images except Figure 5j-l) has the largest spectral value. The spectral value ranges of Band1-3 are narrow, but others are wide, which is highly consistent with the Sentinel-2 image (Figure 5m). Unlike the original GF-2 image whose spectral value peaks in Band2 (Green), Sentinel-2 and all Sentinel-like images have a higher spectral pixel values in the blue band than in the red band, which may be due to the low resolution of Sentinel-2, small crown width, needle-shaped leaves, and large gaps in the Chinese fir woodland. The spectral dispersion of the GS and NND fusion images in the visible light and vegetation red edge bands is better than that of the original Sentinel-2, especially that in Figure 5b,c,f. The spectral distribution of WRM fusion image Figure 5g-i is basically the same as that of Sentinel-2 image. For the BT fusion image Figure 5j-l, the spectral distribution ranges of visible light and vegetation red band are very close. As Figure 6 shows, almost all Sentinellike spectral curves are similar to Sentinel-2, except for the BT fusion images. BT is a simple fusion method that decomposes multi-spectral image pixels into colors and brightness, and then multiplies them with high-resolution images. It can only fuse 3 multi-spectral images at a time. This may be the reason of the distortion in the fusion image Figure 5j-l.
Yang et al. [64] compared the GS, NND and WRM fusion methods, and found that the change trend of the spectral curve of the images fused by three methods are basically the same as that of the original images, but the WRM fusion image retains much more spectral information than other images. This result is consistent with our research.

Selecting Best Fused Image for Forest AGB Estimation
The RF feature selection method is used to quickly screen feature variables from the 14 datasets, including GF-2, Sentinel-2, and the 12 fused Sentinel-like images. By doing this, the feature redundancy and model computation load can be effectively reduced. For each data scenario, that is, 14 groups of feature variables are selected from the above 14 datasets, we use RFR algorithm to establish AGB estimation model. The minimum estimation RMSEr is obtained for image evaluation ( Figure A4).

Selecting Best Fused Image for Forest AGB Estimation
The RF feature selection method is used to quickly screen feature variables fro 14 datasets, including GF-2, Sentinel-2, and the 12 fused Sentinel-like images. By this, the feature redundancy and model computation load can be effectively reduce each data scenario, that is, 14 groups of feature variables are selected from the abo datasets, we use RFR algorithm to establish AGB estimation model. The minimum mation RMSEr is obtained for image evaluation ( Figure A4).
Among the 12 fused Sentinel-like images, NND_B3 (0.2896), GS_B2 (0.3007) GS_B3 (0.3085) have the lowest RMSEr. We compare 14 RMSEr values by the one-sa T test method, and the results show that NND_B3, GS_B2, and GS_B3 images are si cantly different from other images at the 0.05 level, with the p values of 0.000, 0.000 0.021, respectively. The estimation error of these three images is significantly lowe that of other images. The GS method was used to fuse GF-2 multispectral and Lan images in the GSV estimation study of Chinese pine and larch in North China [24 results show that the fusion images obtained based on the B2 and B3 bands of GF-2 im have higher GSV estimation accuracy than other images. In that study, the stepwi gression analysis is used for feature selection, which is different from the RFRmethod, but the results of image data source selection are basically the same, whic ther confirms the feasibility of the improved method. Table 5 shows the normalized statistics of five evaluation indicators of the fuse ages. Image NND_B3 has the values of entropy, standard deviation and average gra significantly different from other images at the 0.05 level, indicating that the imag more information, better texture details and spatial information. Furthermore, its R estimated based on the RF regression algorithm is the lowest, indicating that the i has good forest AGB estimation accuracy. The RMSEr of GS_B2 and GS_B3 are also but are greater than that of NND_B3. Although WRM_B1 has good results in mean, s ard deviation, and entropy, it has the largest RMSEr. Comprehensively considerin information volume, quality and estimation error, we select NND_B3 for the forest modeling and estimation. GS_B2 and GS_B3 images will also be processed to com with NND_B3 images in Section 4.5. Among the 12 fused Sentinel-like images, NND_B3 (0.2896), GS_B2 (0.3007), and GS_B3 (0.3085) have the lowest RMSEr. We compare 14 RMSEr values by the one-sample T test method, and the results show that NND_B3, GS_B2, and GS_B3 images are significantly different from other images at the 0.05 level, with the p values of 0.000, 0.000, and 0.021, respectively. The estimation error of these three images is significantly lower than that of other images. The GS method was used to fuse GF-2 multispectral and Landsat 8 images in the GSV estimation study of Chinese pine and larch in North China [24]. The results show that the fusion images obtained based on the B2 and B3 bands of GF-2 images have higher GSV estimation accuracy than other images. In that study, the stepwise regression analysis is used for feature selection, which is different from the RFR-based method, but the results of image data source selection are basically the same, which further confirms the feasibility of the improved method. Table 5 shows the normalized statistics of five evaluation indicators of the fused images. Image NND_B3 has the values of entropy, standard deviation and average gradient significantly different from other images at the 0.05 level, indicating that the image has more information, better texture details and spatial information. Furthermore, its RMSEr estimated based on the RF regression algorithm is the lowest, indicating that the image has good forest AGB estimation accuracy. The RMSEr of GS_B2 and GS_B3 are also low, but are greater than that of NND_B3. Although WRM_B1 has good results in mean, standard deviation, and entropy, it has the largest RMSEr. Comprehensively considering the information volume, quality and estimation error, we select NND_B3 for the forest AGB modeling and estimation. GS_B2 and GS_B3 images will also be processed to compare with NND_B3 images in Section 4.5.

Selection of Optimal Feature Combination from the Fused Image
The KNN-based and RFR-based feature combination optimization methods are used to select the optimal feature variable combination from the three feature datasets (F1, F2, and F3). Six feature combinations of image NND_B3 are shown in Table 6. For feature set F1, the KNN-based and the RFR-based methods select four and three feature variables, respectively. Three of the seven feature variables are related to SWIR2, indicating that the short-wave infrared band is sensitive to forest vegetation. In F2 and F3, most of the features selected by the two methods are texture feature variables, indicating that the texture feature has a very close relationship with the forest AGB. Table 5. Statistics of normalization for fusion image evaluation index. The RMSEr is calculated by the method in Section 3.3.2. The normalization formula is X' = (X − min)/(max − min), where X' is the normalized data, X is the original data, and max and min are the maximum and minimum values of the original dataset, respectively. The data marked with the symbol * indicates significant difference from other data at the 0.05 level, indicating that these marked images are superior to other unmarked images in corresponding evaluation indexes.

The AGB Estimation Result Analysis
The AGB estimation results of the 4 regression algorithms using the optimal feature combination are shown in Table A1. The KNN-based and RFR-based methods are used in all data scenarios, but only the better results are shown in Table A1. For comparison, three data sources, GF-2, Sentinel-2, and NND_B3, are used to estimate the forest AGB.
The stacking algorithm has good performance in all data scenarios. Take NND_B3 as an example. In the F2 feature set, the RMSEr of the stacking algorithm is 22.09%, which is lower than SVR (24.96%), KNN (24.55%), and RF (22.21%). In the F3 feature set, the RMSEr of stacking (15.53%) is lower than SVR (19.88%), KNN (17.11%), and RF (19.78%) at the statistical level of 0.05. Similarly, in the GF-2 and Sentinel-2 data scenarios, the RMSEr of the stacking algorithm is lower than the other three algorithms, and R 2 , adjusted R 2 , and MAE also have relatively better results. This result is consistent with the performance of the stacking algorithm in the reference [24]. The stacking algorithm used in this study comprehensively considers the coupling and complementarity between models, and optimizes the regression model through the iterative selection of models and the automatic adjustment of hyperparameters, so it can significantly improve the accuracy and stability of AGB estimation of Chinese fir plantations.
Luo et al. [52] studied AGB estimation of the forests in Northeast China based on the Ninth National Forest Continuous Inventory data and Landsat OLI images. They used the recursive feature elimination for feature selection and the categorical boosting as the regression algorithm, and achieved the highest accuracy for coniferous forest, with the RMSE of 26.54 Mg/ha. Coniferous forests in northern China are mainly pine plantations, which are biologically similar to the Chinese fir plantations in southern China. The pine plantations in north China usually have smaller per unit AGB values than Chinese fir plantations in south China [24]. In addition, the topography of the planted forests in northern China is flatter. Therefore, theoretically, the accuracy of AGB estimation of the planted coniferous forests in northern China is usually higher than that in southern China. However, the lowest RMSE and RMSEr value (26.54 Mg/ha, 25.62%) achieved by Luo et al. [52] is larger than our results (15.79 t/ha, 15.53%). On the one hand, their method ignores the ensemble of regression models and the combined effect of feature variables, though they used a large number of field survey plots. On the other hand, they only used Landsat images. Zhang et al. evaluated eight machine learning regression algorithms for estimating forest AGB using satellite remote sensing data and multiple auxiliary data [54]. They concluded that the categorical boosting (CatBoost) algorithm has the best accuracy among these eight algorithms, with the R 2 (0.72), RMSE (45.63 Mg/ha), and RMSEr (25%). However, the CatBoost achieved poor performances in estimating the AGB of evergreen needleleaf forests, as the RMSEr is larger than 60%. The poor accuracy could be resulted from the underestimation of evergreen needleleaf forest samples due to the saturation problems.

The AGB Estimation Ability of Different Image Data Source
The NND_B3 obtained by fusing the red band of GF-2 with Sentinel-2 by the NND method has the advantages of both GF-2 and Sentinel-2. GF-2 and Sentinel-2 have similar AGB estimation performance in F1 feature sets. However, for feature sets F2 and F3, the RMSEr of NND_B3, GF-2 and Sentinel-2 are different (Table A1). For the F2 feature set, the optimal RMSEr of Sentinel-2 is 0.0288 larger and R 2 is 0.1704 smaller than that of NND_B3. For the F3 feature set, the optimal RMSEr of NND_B3 is lower than Sentinel-2 by 0.0515, and R 2 is larger by 0.2343. Using NND_B3, the estimation accuracy is greatly improved. As the scatter plot of Figure A5 shows, NND_B3 has more concentrated distribution, clearer fitting trend, and higher correlation coefficient r than other data. We also compared the estimation accuracy of using GS_B2 (RMSEr, 0.1804) and GS_B3 (RMSEr, 0.1893) and found their accuracy is lower than using NND_B3(RMSEr, 0.1553) (Figure 7).

The Best Feature Selection Method for Different Data Scenarios and Different Estimation Models
Different feature selection methods are applicable to different estimation models. In this study, KNN-based and RFR-based were used to select the optimal combination of feature variables for F1, F2, and F3 feature sets. The estimation results of 36 models indicate that (Table A1), for the F2 feature set of Sentinel-2 and NND_B3, the RFR-based feature variable selection method is better than the KNN-based in all models. In addition, as Table A1 and Figure A4 show, RFR-based has a lower RMSEr value in the F1 feature set of GF-2, Sentinel-2, and NND_B3 images than the RF importance index method. This means that the RFR-based method has better performance than traditional RF-base method in the feature selection of forest AGB estimation. For most data scenarios, when the RF regression algorithm is used, the RFR-based performs better than the KNN-based. The KNN-based is better than the RFR-based for the SVR, KNN, and stacking models of the F1 and F3 feature variable set. This result indicates that the RFR-based method can be combined with the KNN-based method in different model application scenarios.
Remote Sens. 2021, 13, x FOR PEER REVIEW 1 Figure 7. The estimation results of three images by the stacking algorithm.

The Best Feature Selection Method for Different Data Scenarios and Different Estimati Models
Different feature selection methods are applicable to different estimation mode this study, KNN-based and RFR-based were used to select the optimal combinati feature variables for F1, F2, and F3 feature sets. The estimation results of 36 models cate that (Table A1), for the F2 feature set of Sentinel-2 and NND_B3, the RFR-base ture variable selection method is better than the KNN-based in all models. In additi Table A1 and Figure A4 show, RFR-based has a lower RMSEr value in the F1 featu of GF-2, Sentinel-2, and NND_B3 images than the RF importance index method means that the RFR-based method has better performance than traditional RF method in the feature selection of forest AGB estimation. For most data scenarios, the RF regression algorithm is used, the RFR-based performs better than the KNN-b The KNN-based is better than the RFR-based for the SVR, KNN, and stacking mod the F1 and F3 feature variable set. This result indicates that the RFR-based method c combined with the KNN-based method in different model application scenarios.

AGB Estimation Performance of Different Feature Sets
The vegetation index and texture feature variables of optical remote sensing im can be used for forest classification and structural parameter prediction [38]. The AGB estimation performance of different feature sets varies significantly (Table A performs the best, F2 is the second, and F1 is the worst. Taking the NND_B3 image example, when the RF and stacking estimation models are used, the RMSEr of F3 is 0 and 0.0656 less than that of F2, respectively, and the MAE is 1.78 t/ha and 5.0 t/ha respectively. Similarly, compared with the RESEr of F1, the RMSEr of F2 is 0.022 0.0269 less, respectively, and the MAE is 1.83 t/ha and 3.78 t/ha smaller, respectively. parison shows that the estimated AGB results of F3 are more correlated with the mea AGB values ( Figure A5). The average value of r is 0.6868, and the maximum valu reach 0.8427. F1 performs poorly, with the mean value of r being 0.4143 and the mini value being 0.1614. Figure A6 shows the variation range and trend of AGB estimation deviation o ferent feature sets of NND_B3 image. The overall estimation deviation of F3 is rela low. Larger RMSErs are mainly of the samples whose measured AGB exceeds 140 and there are only 11 samples with the absolute deviation exceeding 20 t/ha. F2 a

AGB Estimation Performance of Different Feature Sets
The vegetation index and texture feature variables of optical remote sensing images can be used for forest classification and structural parameter prediction [38]. The forest AGB estimation performance of different feature sets varies significantly (Table A1). F3 performs the best, F2 is the second, and F1 is the worst. Taking the NND_B3 image as an example, when the RF and stacking estimation models are used, the RMSEr of F3 is 0.0243 and 0.0656 less than that of F2, respectively, and the MAE is 1.78 t/ha and 5.0 t/ha less, respectively. Similarly, compared with the RESEr of F1, the RMSEr of F2 is 0.0222 and 0.0269 less, respectively, and the MAE is 1.83 t/ha and 3.78 t/ha smaller, respectively. Comparison shows that the estimated AGB results of F3 are more correlated with the measured AGB values ( Figure A5). The average value of r is 0.6868, and the maximum value can reach 0.8427. F1 performs poorly, with the mean value of r being 0.4143 and the minimum value being 0.1614. Figure A6 shows the variation range and trend of AGB estimation deviation of different feature sets of NND_B3 image. The overall estimation deviation of F3 is relatively low. Larger RMSErs are mainly of the samples whose measured AGB exceeds 140 t/ha, and there are only 11 samples with the absolute deviation exceeding 20 t/ha. F2 and F1 have 19 and 24 samples with the absolute deviation exceeding 20 t/ha, respectively, and 5 samples of F1 have the absolute deviation larger than 40 t/ha. Figure 8(a1-c3) are the AGB spatial distribution map of Chinese fir plantation estimated by the Stacking algorithm using the GF-2, Sentinel-2 and NND_B3, respectively. Figure 8(a1,b1,c1) are the estimation results of the F1 feature set of GF-2, Sentinel-2 and NND_B3, respectively, whose AGB values are between 90 t/ha and 110 t/ha, generally low. This indicates that the estimated AGB value is less saturated and overestimated for low values. The results of the F2 feature set (Figure 8(a2,b2,c2)) have less overestimation, and the estimated values range between 69 t/ha to 158 t/ha. The results of the F3 feature set are the best, with AGB ranging between 45 t/ha and 176 t/ha. The green area of Figure 8(c3) is larger than Figure 8(a3,b3), indicating more AGB estimates are between 135t/ha and 176t/ha. In short, NND_B3 has better estimation results than GF-2 and Sentinel-2, and it improves the saturation significantly. The F3 feature set with multi-window texture factors is better than F1 and F2 in the lower AGB area. It has better performance, and supresses the overestimation, indicating that the multi-window texture factor can suppress the overestimation in the low-value area of AGB. Figure 8(d1-d3) show the AGB result distribution map estimated by SVR, KNN, and RF model using the F3 feature set of NND_B3 image. The obtained AGB values range between 44 t/ha and 168 t/ha, 49 t/ha and 176 t/ha, and 59 t/ha and 161 t/ha, respectively. The SVR model has higher accuracy in the low-value area, but lower accuracy in the high-value (green) area than the KNN and RF models. In the southern area, Figure 8(c3,d2) and d3 all have higher AGB estimates, but Figure 8(d1) has low AGB, which also indicates that the SVR model has the underestimation problem in the area with higher AGB values. Finally, compared with the estimation results of Figure 8(c3,d1-d3), Figure 8(c3) has more high value areas and wider AGB value range, indicating that the stacking model has stronger generalization ability and higher accuracy than SVR, KNN, and RF. Zhao et al. [65] studied the AGB estimation data saturation problem using Landsat TM (Thematic Mapper) images for different vegetation types and obtained the forest biomass saturation values of 159 Mg/ha for pine (Pinus Massoniana) plantations forests in Eastern China. Gao et al. [19] compared several models for AGB estimation in subtropical forests, and found that the RF algorithm is not suitable for AGB prediction when the AGB values are too small (<40 Mg/ha) or too large (>160 Mg/ha). This is consistent with the AGB estimation result based on the RFR model in this study.

Limitations and Future Works
This study has proved the superiority of multi-spectral fusion image combined with stacking integrated modeling method in estimating the AGB of Chinese fir plantation, but there are some limitations in the application. First of all, optical remote sensing images are usually polluted by clouds. This affects the imaging quality of remote sensing images to a certain extent, and even significantly reduces the availability of image data. Therefore, the difference in imaging time for GF-2 and Sentinel-2 brings certain uncertainty to image fusion processing. Secondly, as shown in Table 2, more than half of the Chinese fir plantation sample plots are not really mature, meaning that tree growth that occurred between 2016 and 2017 has not bet accounted for. Third, the stacking integration algorithm usually has a great computation load [24,55], which should be reduced. The categorical boosting regression algorithm and extreme gradient boosting regression algorithm works well in classification and regression prediction [48,52,53], so combining them with the stacking ensemble algorithm may improve the forest AGB estimation. Remote Sens. 2021, 13, x FOR PEER REVIEW 20 of 30

Conclusions
In this study, the RF-S method was proposed for estimating the AGB of Chinese fir plantations in south China. The results demonstrate the superiority of fusing GF-2 multispectral and Sentinel-2 data, as well as the potential of the improved feature combination optimization method at regional scales. The stacking generalization method based on the fused image (NND_B3) has higher estimation accuracy and saturation than other images and regression models. The achieved adjusted R 2 and RMSEr are 0.6306 and 15.53%, respectively. This study uses both vegetation index and texture feature factors of multiple window sizes as the input feature variables for the training model, which can provide higher accuracy and data saturation for Chinese fir plantation AGB mapping. However, the proposed RF-S strategy has huge computation load, due to forest plot data collection, multi-spectral image fusion processing, feature variable combination optimization, and integration of multiple models. Thus, methods for improving the computation efficiency should be developed before this strategy is widely applied. This research provides some insights for forest AGB estimation research based on remote sensing images and sample plots data modeling.

Data Availability Statement:
The observed GSV data from the sample plots and spatial distribution data of forest resources presented in this study are available on request from the corresponding author. Those data are not publicly available due to privacy and confidentiality. The Sentinel-2 images of the L1-level product were obtained from Copernicus data center website at https://scihub. copernicus.eu/ (accessed on 17 March 2020). The GF-2 images are available from China Centre for Resources Satellite Data and Application website at http://www.cresda.com/CN/ (accessed on 10 May 2020). The DEM data were downloaded from geospatial data cloud (http://www.gscloud.cn/ (accessed on 15 June 2020)).

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

Appendix A
Supplementary notes: Because BT fusion method can only fuse 3 bands at a time, in this study, based on BT fusion method, we only get 9 fused bands. The regression modeling algorithms used in this study, including RF, KNN, SVR, and stacking, are implemented by calling the machine learning toolkit based on Python 3.7 programming language. In addition, the image GLCM calculation was performed based on the "Co-occurrence Measures" function of Envi 5.3 software, and image sharpening fusion was performed using Envi 5.3 and Erdas 9.2 software.
Supplementary notes: Because BT fusion method can only fuse 3 bands at a time, in this study, based on BT fusion method, we only get 9 fused bands. The regression modeling algorithms used in this study, including RF, KNN, SVR, and stacking, are implemented by calling the machine learning toolkit based on Python 3.7 programming language. In addition, the image GLCM calculation was performed based on the "Co-occurrence Measures" function of Envi 5.3 software, and image sharpening fusion was performed using Envi 5.3 and Erdas 9.2 software.   . Figure A2. The RFR-based method for optimizing the feature variable combination. (a) Selection of the best fused image; (b) Feature variable selection based on combinatorial optimization method; (c) Basic flow of random forest regression algorithm.    Figure A5. The scatter graphs between the observed and estimated AGB values of the Chinese fir plots using three image datasets and four estimation models: (a-d), (e-h) and (i-l) are the AGB estimated by the GF-2, Sentinel-2, and NND_B3 image using the SVR, KNN, RF, and stacking algorithm, respectively. Blue, orange, and gray represent F1, F2, and F3 feature variable sets, respectively. Figure A5. The scatter graphs between the observed and estimated AGB values of the Chinese fir plots using three image datasets and four estimation models: (a-l) are the AGB estimated by the GF-2, Sentinel-2, and NND_B3 image using the SVR, KNN, RF, and stacking algorithm, respectively. Blue, orange, and gray represent F1, F2, and F3 feature variable sets, respectively. Remote Sens. 2021, 13, x FOR PEER REVIEW 26 of 30 Figure A6. The estimation bias of the Stacking algorithm based on the NND_B3 image :(a-c) are F3, F2, and F1 feature variable sets, respectively. The red dotted line is the deviation change trend line. Figure A6. The estimation bias of the Stacking algorithm based on the NND_B3 image: (a-c) are F3, F2, and F1 feature variable sets, respectively. The red dotted line is the deviation change trend line. Table A1. The summary of assessment indicators for AGB estimation results of Chinese fir plantations based on nine image feature data scenarios and four estimation algorithms. The estimation results are the better ones of the feature screened by the KNN-based and RFR-based algorithms. The data scene with light green shade in the table indicates that the estimation result of feature variables selected by the RFR-based algorithm is better than the KNN-based method. However, the KNN-based algorithm is better in other data scenarios.