Co-Kriging-Guided Interpolation for Mapping Forest Aboveground Biomass by Integrating Global Ecosystem Dynamics Investigation and Sentinel-2 Data

Abstract

Mapping wall-to-wall forest aboveground biomass (AGB) at large scales is critical for understanding global climate change and the carbon cycle. In previous studies, a regression-based method was commonly used to map the spatially continuous distribution of forest AGB with the aid of optical images, which may suffer from the saturation effect. The Global Ecosystem Dynamics Investigation (GEDI) can collect forest vertical structure information with high precision on a global scale. In this study, we proposed a collaborative kriging (co-kriging) interpolation-based method for mapping spatially continuous forest AGB by integrating GEDI and Sentinel-2 data. First, by fusing spectral features from Sentinel-2 images with vertical structure features from GEDI, the optimal estimation model for footprint-level AGB was determined by comparing different machine-learning algorithms. Second, footprint-level predicted AGB was used as the main variable, with rh95 and B12 as covariates, to build a co-kriging guided interpolation model. Finally, the interpolation model was employed to map wall-to-wall forest AGB. The results showed the following: (1) For footprint-level AGB, CatBoost achieved the highest accuracy by fusing features from GEDI and Sentinel-2 data (${\mathrm{R}}^2$ = 0.87, RMSE = 49.56 Mg/ha, rRMSE = 27.06%). (2) The mapping results based on the interpolation method exhibited relatively high accuracy and mitigated the saturation effect in areas with higher forest AGB (${\mathrm{R}}^2$ = 0.69, RMSE = 81.56 Mg/ha, rRMSE = 40.98%, bias = −3.236 Mg/ha). The mapping result demonstrates that the proposed method based on interpolation combined with multi-source data can be a promising solution for monitoring spatially continuous forest AGB.

1. Introduction

Forests are an important part of terrestrial ecosystems and play a significant role in the global carbon budget [1,2]. Forest aboveground biomass (AGB) is an important indicator for judging the carbon resources and sequestration capacity of terrestrial ecosystems [3,4]. Accurate measurements of forest AGB at large scales are of great significance for quantifying terrestrial carbon stocks and monitoring forest ecosystem productivity [5,6]. However, it is not easy to map high-precision wall-to-wall forest AGB, and the generation of large-scale forest AGB products is still a current research hotspot [7,8,9].

Traditional forest field surveys are labor-intensive and time-consuming, so it is essential to utilize available remote-sensing data for the estimation of large-scale forest AGB. Remote-sensing data used for AGB estimation mainly include multispectral remote sensing (MRS), synthetic aperture radar (SAR), and light detection and ranging (LiDAR) [10,11,12]. Spectral reflectance information of the vegetation surface can be collected by optical images, which enables quantitative inversion of forest structural parameters [13]. Fardin et al. [14] modeled AGB in sparse Mediterranean forest ecosystems using band information and spectrally derived vegetation indices from Sentinel-2 and then found that the infrared percentage vegetation index (IPVI) and normalized difference vegetation index (NDVI) had the highest correlation with AGB. Unfortunately, optical sensors may have strong saturation effect in dense forests [15,16], and their estimated forest AGB products are usually fraught with uncertainty issues. SAR sensors have a greater penetration capability than optical sensors. The information below the vegetation canopy can be obtained [17], and thus, it is more suitable for the estimating of forest AGB. Zhang et al. [18] extracted the backscattering coefficients of four polarization modes (VV, VH, HV, and HH) from Sentinel-1 to estimate the AGB and found that the VV polarization mode had the highest accuracy. However, the backscattering characteristics of SAR data are easily affected in regions with complex terrain, potentially impacting the accuracy of AGB estimation [19]. LiDAR, as an active remote-sensing technology, can obtain accurate forest structural information by using a focused short-wavelength laser pulse [20,21]. This technology plays a significant role in improving the accuracy of forest AGB prediction. For example, Ihor et al. [22] improved the estimation of AGB of Pinus sylvestris in urban forest by using airborne LiDAR and digital hemispherical photography (DHP) techniques, with an ${\mathrm{R}}^2$ of 0.98 and an RMSE of 2.97 t/ha. Unmanned aerial vehicle LiDAR (UAV-LiDAR) was also used in conjunction with Sentinel-2, substantially improving the accuracy of bamboo forest AGB estimation [23]. These study results indicated that a relatively higher accuracy in AGB estimation can be achieved by utilizing LiDAR technology. However, it is impractical to use the airborne platform to collect large-scale forest structure observations considering the high cost of acquiring LiDAR data [24,25]. To obtain vertical structure information of forests in large-scale research, studies based on spaceborne LiDAR data have begun to emerge [26,27]. Petri et al. [28] proposed a hierarchical hybrid inference approach for uncertainty quantification of the average aboveground biomass density (AGBD) estimated directly from a sample of spaceborne LiDAR profiles. The results show that the reference estimate was within the 95% confidence interval of the hierarchical hybrid estimate, which provided an effective way to quantify the uncertainty of estimated AGBD. However, all spaceborne LiDAR systems share a common problem of having discrete ground sampling footprints, and they cannot directly provide wall-to-wall forest biomass surveys [29].

To overcome the spatial discontinuity issue of spaceborne LiDAR, researchers have attempted to combine spaceborne LiDAR with optical imagery to map wall-to-wall forest AGB products [30,31]. In these studies, vertical structure information from spaceborne LiDAR and spectral reflectance information are used as predictors for building regression models. For example, ICESat-2/ATLAS (Ice, Cloud, and Land Elevation Satellite-2, Advanced Topographic Laser Altimeter System) and Landsat 8 OLI data were employed to estimate forest AGB in the southeastern region of Texas, with an ${\mathrm{R}}^2$ of 0.58 and an RMSE of 23.89 Mg/ha [32]. For different forest types, GEDI, Sentinel-1, and Sentinel-2 were employed to estimate and map the AGB, demonstrating great potential for regional-scale AGB mapping [33]. Wang et al. [34] proposed a random-forest-based downscaling method to map forest height and biomass at a 15 m resolution by integrating Landsat 8 OLI and ICESat-2 LiDAR data. Subrata et al. [35] mapped forest canopy height by integrating ICESat-2 and Sentinel-1 data and used canopy height information along with spectral variables derived from Sentinel-2 data to estimate forest AGB in the foothills of the Himalayas in north-west India. For these studies, spectral features from optical images are usually fed into a machine-learning model to map wall-to-wall forest AGB. As such, these biomass products are still affected by saturation effects. Furthermore, the spatial and temporal heterogeneity of AGB is generally ignored in these studies [36].

Geostatistical methods can explain spatial heterogeneity and correlations among variables, such as biomass, canopy height, and canopy cover [37]. Its spatial interpolation technique uses observations from known locations to predict values at unknown locations [38], further building spatial prediction models that have been widely applied in various fields [39,40,41]. Spaceborne LiDAR data with a wide distribution can obtain relatively high-precision vertical structure information of vegetation, facilitating the application of interpolation methods [42]. However, forest biomass interpolated from the spaceborne LiDAR platform may display a strip effect, as spaceborne LiDAR footprints are generally evenly distributed along ground tracks with an interval [43]. For geostatistical co-kriging, both the main variable and covariates are jointly used to capture spatial heterogeneity, which can mitigate the strip effect and reduce estimation error. Therefore, by utilizing spaceborne LiDAR data with densely distributed footprints, the interpolation method may be a promising solution for mapping spatially continuous forest AGB.

In this study, we propose a method for mapping high-quality and spatially continuous forest AGB by integrating multi-source data with collaborative kriging (co-kriging) interpolation. First, parametric and non-parametric methods are utilized for feature selection of GEDI and Sentinel-2 data, providing reliable independent variables for AGB modeling. Second, CatBoost, random forest (RF), and multiple linear regression (MLR) are employed to determine the optimal forest AGB estimation model at the footprint level. Then, the optimal model is selected to predict AGB within all footprints in the study area. Finally, based on the selected semivariogram, a co-kriging spatial prediction model is established for mapping wall-to-wall forest AGB. This study mitigated the saturation effect in areas with higher AGB to some extent and more accurately mapped the spatially continuous forest AGB at the county level.

2. Materials

2.1. Study Area

The study area is in Sonoma County, California, United States, situated between the Pacific Ocean and Napa Valley (121.76°W–122.46°W, 38.11°N–38.85°N). Sonoma County is located in the north of San Francisco with a total area of about 4580 km², of which 89.12% is land and 10.88% is water. The area’s climate is characterized as Mediterranean, with annual precipitation averaging approximately 500 mm in the southeastern county and ranging from 700 to 1000 mm in the central and northern valley areas. Meanwhile, the average monthly temperature ranges from 7.3 to 22.6 °C throughout the year [44]. The topography of Sonoma County is complex and mainly consists of several mountains and hills, with elevations ranging from −20 to 1366 m (see Figure 1).

Vegetation in the study area is diverse and primarily includes coniferous forests, broadleaf forests, mixed forests, and shrubs. Among these, oak (broadleaf forests) is the dominant forest type, making up 36% of the forested area. Coniferous forests are dominated by firs, pines, and sequoias, which account for over 38%. Mixed forests and shrubs follow, exhibiting a large gradient of overall forest structure [45,46].

2.2. Airborne Laser Scanning (ALS)

The NASA Carbon Monitoring System (CMS) project provided the AGB dataset for Sonoma County, which is free and publicly available (https://daac.ornl.gov/ (accessed on 10 August 2023)). The dataset is provided in GeoTIFF format with a spatial resolution of 30 m. Forest AGB was generated using a combination of airborne LiDAR data and field plot measurements with a parametric modeling approach. The field estimates of AGB from 194 sample plots were related to height and cover metrics derived from ALS, and the uncertainty of the parametric model used to generate the AGB map was also quantified, resulting in a relatively high-accuracy AGB map [47].

The point density of the LiDAR campaign was approximately 14 points per square meter across 44 flight lines, covering an area of 440,000 ha. The dataset was used as the reference AGB in this study for constructing and evaluating AGB models at the footprint level, as well as for verifying the accuracy of AGB mapping.

2.3. GEDI

GEDI data were downloaded from the earth data platform of the National Aeronautics and Space Administration (NASA) (https://search.earthdata.nasa.gov/ (accessed on 15 August 2023)). The total width of GEDI laser tracks is 4.2 km, with an along-track spacing of 60 m, a cross-track spacing of 600 m, and an average footprint size of 25 m. L2A and L2B data from GEDI used in this study are canopy structure products at the footprint level, with a spatial resolution of 25 m. Specifically, the L2A product contains information derived from the geolocated GEDI return waveforms, including ground elevation, highest and lowest surface return elevations, and relative height metrics, while biophysical information is derived from the L2B product [48,49]. From the L2 products, the vertical distribution information of the vegetation canopy can be easily extracted and used for the estimation of forest AGB. The detailed description of the L2 product used in this study is shown in Table 1.

2.4. Sentinel-2

Sentinel-2 is a high-resolution multispectral imaging satellite with ground resolutions of 10, 20, and 60 m [50]. The satellite covers 13 spectral bands ranging from visible to near-infrared to short-wave infrared. Currently, Sentinel-2A is the only optical satellite that includes three red-edge bands, which are effective for monitoring vegetation health and growth [51]. Therefore, the spectral information of vegetation from Sentinel-2A was utilized in this study to construct forest AGB estimation models.

2.5. Ancillary Data

The Shuttle Radar Topography Mission (SRTM) globally acquired high-quality elevation data, providing digital elevation datasets with a spatial resolution of 30 m [52]. In this study, we utilized the SRTM dataset, accessed through the Google Earth Engine (GEE) platform, to extract ground elevation, slope, and aspect, which can effectively reduce the topography’s impact on biomass estimation [53].

The Global 30 m Fine Ground Cover 2020 product [54] can meet application analysis needs on a global scale. The product, known for its high accuracy and detailed classification of surface vegetation types, was used in this study to classify the vegetation types on the surface of the study area.

3. Methodology

In this study, we propose a method for mapping high-precision, wall-to-wall forest AGB at the county level based on co-kriging interpolation of GEDI and Sentinel-2A data. Figure 2 presents an overview of the AGB mapping process; it can be divided into the following three steps.

The first step is data preprocessing and feature extraction. Low-quality GEDI footprints and Sentinel-2 images are removed by data preprocessing. Then, three types of footprint-level features are extracted from different data sources based on the positions of the GEDI footprints. The second step is to build AGB regression models at the footprint level. Different screening methods are employed to filter the input variables for parametric and non-parametric estimation models. The estimation models are established using the CatBoost, RF, and MLR algorithms, and ${\mathrm{R}}^2$, RMSE (Mg/ha), and rRMSE (%) are compared to determine the optimal model. The third step involves mapping forest AGB. A co-kriging interpolation model is built and then used to map the wall-to-wall forest AGB across the study area. During the interpolation process, the sample points of biomass can be constructed through the generated footprint-level AGB.

3.1. Data Preprocessing

The data preprocessing mainly involves the removal of invalid GEDI footprints and the preprocessing of Sentinel-2 images. Due to the noise influence on GEDI data, the estimation of AGB would be affected if the data are used directly. High-quality GEDI footprints were obtained by filtering out invalid ones based on the experience of previous studies [42,55,56]. As a result, a total of 77,310 valid footprints were identified in the study area. The valid footprints filtered for urban and vegetation areas are shown in Figure 3. The specific filtering conditions applied in this study are as follows:

• Lon_lowestmode, Lat_lowestmode, and shot_number: These values indicate the longitude, latitude, and identification number of the footprints, which can be used to find the footprints’ locations.

• quality_flag: The value indicates the quality of the footprint. If the quality_flag is 1, it means that the footprint meets the criteria based on energy, sensitivity, amplitude, and real-time surface tracking quality, and thus is retained.

• degrade_flag: If the value is 1, it means that the state of the pointing or geolocated information is degraded; thus, only the footprints with degrade_flag = 0 are retained.

• Sensitivity: The probability that the echo signal reaches the ground from the top of the canopy, with a value greater than or equal to 0.9 representing good spot quality. The sensitivity thresholds varied for different land cover types, with a threshold of ≥ 0.95 set in this study.

Furthermore, Sentinel-2 images were processed using the GEE platform. To enhance image quality and thereby improve the AGB inversion accuracy, only images with less than 5% cloud cover were selected based on the cloud removal algorithm. Since the different bands of Sentinel-2 have varying resolutions, the bands used in this study were resampled to a uniform resolution of 30 m before calculating the vegetation index.

3.2. Feature Variable Extraction

In this study, the features used for constructing AGB models were primarily divided into three categories: GEDI features, Sentinel-2 image features, and SRTM DEM terrain factors. The total number of these features is 230 (see Table 2).

GEDI features refer to the forest vertical structure parameters of the L2A and L2B products, with a total number of 39. L2A is primarily a product related to tree height, which indicates the maturity of forests. Forests with taller trees are usually more mature, and thus have relatively higher biomass. The canopy cover and vertical profile parameters from L2B serve as biophysical metrics, which have a significant impact on the forest AGB prediction (see Table S1).

Different spectral bands have varying reflectivity characteristics for different objects. Therefore, 10 bands of the Sentinel-2 images (B2, B3, B4, B5, B6, B7, B8, B8a, B11, and B12) were first selected for forest AGB prediction. The vegetation index is one of the most important variables for estimating forest AGB because it provides information on vegetation growth status, chlorophyll content, and vegetation cover (see Table S2) [57]. Texture features, which are the variations in the grayscale of images, reflect the surface characteristics of ground objects and are crucial for image interpretation. The gray-level co-occurrence matrix (GLCM) is a commonly used method for describing texture by analyzing the spatial correlation of grayscale values, which can capture the detailed features within the forest. In this study, a total of 188 image features were extracted for forest AGB prediction, including 10 spectral bands, 8 vegetation indices, and 170 texture features (see Table S3).

Terrain factors play an important role in estimating forest biomass. In this study, three terrain factors were used: elevation, slope, and aspect (see Table S4). Elevation can influence climate and soil conditions, and therefore the growth of forests. Similarly, slope and aspect can affect rates of soil erosion and water distribution, thus influencing forest growth and biomass accumulation [58].

3.3. AGB Estimation at the Footprint Level

3.3.1. Feature Variable Selection

Too many features increase the complexity of the estimation model, resulting in the risk of overfitting, increased training time, and computational costs [59]. Therefore, it is necessary to perform feature selection before establishing the model. By this means, the efficiency and generalization ability of the model can be effectively improved. The process of feature selection in this study includes two steps: correlation analysis and feature optimization.

Correlation analysis, a statistical method used to measure the degree of relationship between two or more variables, is initially performed in this study by using SPSS 22.0 software to identify features significantly correlated with AGB. Subsequently, two methods are employed to further optimize the independent variables of the prediction model. Specifically, the stepwise method and random forest are used to selectively filter variables in the parametric regression model and the non-parametric regression model, respectively. When using the stepwise method for feature optimization in the MLR model, variables are added one by one. Meanwhile, the variance inflation factor (VIF) is introduced, and the corresponding Akaike Information Criterion (AIC) values are calculated to determine the combination of variables that minimizes the AIC. Random forest calculates the relative importance of each variable for the model’s predictive performance by integrating multiple decision trees, which helps select independent variables with high importance scores. Finally, the optimal combination of features can be derived for further AGB modeling.

3.3.2. Footprint-Level AGB Modeling

To compare the estimation accuracy of different models, footprint-level forest AGB estimation models were established by using three algorithms: CatBoost, RF, and MLR, which were implemented in Python 3.10.

CatBoost is an advanced machine-learning algorithm that employs symmetric decision trees as base learners to fit the residuals [60]. Based on the residuals from the previous tree, each new tree is iteratively optimized. Finally, the predictions of all the decision trees are accumulated to produce the model’s final prediction. The algorithm optimizes issues related to gradient bias and prediction shift and is characterized by fewer parameters, fast convergence, and strong robustness, which enhance the model’s generalization ability [61].

RF is an ensemble learning algorithm based on decision trees [62,63]. For regression tasks, RF employs the bootstrap sampling method to randomly select a subset of samples with replacement from the dataset and constructs decision trees based on a randomly selected subset of features. The results of all the decision trees are weighted and averaged to generate the final prediction for the sample [64]. The algorithm is characterized by fast computational efficiency, good fitting performance, and strong stability, all of which can effectively improve estimation accuracy in modeling [65].

MLR, as a parametric regression model, aims to fit the linear relationship between one dependent variable and two or more independent variables [42]. The general form of MLR is as follows:

(1)$y={\beta }_0+{\beta }_1{x}_1+{\beta }_2{x}_2+{\beta }_3{x}_3+\dots +{\beta }_i{x}_i+\epsilon$

where y is the dependent variable, x is the independent variable, ${\beta }_0$ is the regression constant, $\epsilon$ is the random error, and ${\beta }_1,{\beta }_2,{\beta }_3,{\beta }_i$ are the regression coefficients.

The optimal AGB prediction model is determined by comparing the testing accuracy, modeling efficiency, and complexity of parameter settings of each algorithm, enabling the estimation of forest biomass for all GEDI footprints.

3.4. Mapping Wall-to-Wall Forest AGB Based on the Geostatistical Co-Kriging Method

Geostatistics is commonly used for analyzing spatial data, with the core concept being to predict data at unknown locations by quantifying the spatial relationships (e.g., distance and direction) between sample points [66]. Co-kriging is an advanced spatial interpolation method in geostatistics that can transform a limited number of points into the surface of a region. Therefore, the co-kriging method was employed in this study to map wall-to-wall forest AGB. Firstly, the footprint-level predicted AGB was used as the main variable for co-kriging interpolation. Since rh95 and B12 demonstrated high importance scores in the feature importance analysis and were strongly correlated with AGB, while considering the highly variable topography of the area, these features (rh95, B12, and slope) were selected as covariates. Secondly, the optimal semivariogram was determined through semivariance function analysis to accurately model the spatial correlation among AGB sample points. Finally, the interpolation model was used to generate the wall-to-wall forest AGB.

3.4.1. Semivariance Function

The semivariance function is an important tool in co-kriging, used to quantify and model spatial correlations among sample points. During the interpolation process, these correlations accurately reflect the spatial structure of the data and play a significant role in the allocation of weights, providing more precise predictions for unknown points. To improve the accuracy and effectiveness of interpolation results, it is necessary to fit an empirical semivariogram model (i.e., the visual representation of the semivariance function) to analyze the correlations among sample points and help obtain the appropriate weights for predictions [67]. The common functions for modeling the empirical semivariogram include linear, spherical, exponential, and Gaussian, which are described using the following parameters: nugget, sill, and range [68].

The semivariance function represents the variance of the regionalized variable (i.e., footprint-level AGB in this study) over a certain distance, indicating how the spatial correlation of the variable changes within the range. The formula is as follows:

(2)$\gamma (h)={\displaystyle \frac{1}{2N(h)}}{\displaystyle \sum _{i=1}^{N(h)}{\left[Z(x_i)-Z(x_i+h)\right]}^2}$

where $\gamma (h)$ is the value of the semivariance function, $N(h)$ is the number of pairs of points with a distance equal to h, and $Z(X_i)$ and $Z(x_i+h)$ are the reference values of AGB at point $x_i$ and point $x_i+h$, respectively.

3.4.2. Co-Kriging Interpolation

Through the analysis of the semivariance function, the spatial correlation between AGB sample points can be determined. The $\gamma (h)$ values indicate how the correlation between samples changes with distance, providing the input parameters for weight calculation. According to these findings, co-kriging is then used to interpolate the forest AGB. Co-kriging is an optimal method for the linear unbiased estimation of regionalized variables within a restricted area [37]. By combining covariates with the main variable, this method provides both the spatial autocorrelation modeling of ordinary kriging and the trend information of biomass, thus improving the precision of spatial data interpolation. The formula is as follows:

(3)${\widehat{z}}_0={\displaystyle \sum _{i=1}^n{\lambda }_i{z}_i}+{\displaystyle \sum _{j=1}^m{b}_j{x}_j}$

where ${\widehat{z}}_0$ is the estimated AGB of the unknown point, $z_i$ is the value of the main variable at point $i$, $x_j$ is the value of the covariate at point $j$, and ${\lambda }_i,b_i$ are the weights to be estimated.

By minimizing the variance of the co-kriging estimation error, the weights ${\lambda }_i$and $b_i$can be determined [40]. To ensure the estimates are unbiased under the second-order smooth condition, assume the following equation:

(4)$\left\{\begin{array}{c}\left[\begin{array}{ccc}{\gamma }_{zz}& {\gamma }_{zx}& 1\\ {\gamma }_{xz}& {\gamma }_{xx}& 0\\ 1& 0& 0\end{array}\right]\left[\begin{array}{c}\lambda \\ b\\ \mu \end{array}\right]=\left[\begin{array}{c}{\gamma }_z\left(y_0\right)\\ {\gamma }_x\left(y_0\right)\\ 1\end{array}\right]\\ \begin{array}{c}{\displaystyle \sum _{i=1}^n{\lambda }_i=1}\\ {\displaystyle \sum _{j=1}^m{b}_j=0}\end{array}\end{array}\right.$

where ${\gamma }_{zz}$ is an n-order matrix composed of semivariance values calculated between sample points of the main variable z, ${\gamma }_{xx}$ is an m-order matrix composed of semivariance values calculated between sample points of the covariate x, ${\gamma }_{zx}$ is a matrix composed of semivariance values calculated for the main variable and the covariate, ${\gamma }_z\left(y_0\right)$ is the semivariance matrix calculated between the samples of the main variable z and the prediction location $y_0$, ${\gamma }_x\left(y_0\right)$ is the semivariance matrix calculated between the samples of the covariate x and the prediction location $y_0$, and $\mu$ represents the Lagrange multiplier.

3.5. Accuracy Evaluation

The accuracy of AGB estimation models was assessed using evaluation metrics, including the coefficient of determination ($R^2$), root mean square error (RMSE), and relative root mean square error (rRMSE). $R^2$reflects the goodness of fit; the closer the value is to 1, the better the regression model is fitted. RMSE reflects the dispersion of the samples and is used to measure the deviation between predicted values and actual values; a smaller RMSE indicates a better model. rRMSE, defined as the ratio of RMSE to the true values, reflects the accuracy of the models; a smaller rRMSE value implies higher accuracy. These metrics are widely employed in forest AGB modeling to assess the differences between measured and predicted values [23,68,69]. The calculation formulas are as follows:

(5)$R^2=1-{\displaystyle \frac{{\displaystyle {\sum }_{i=1}^n(y_i-{\widehat{y}}_i)}}{{\displaystyle {\sum }_{i=1}^n(y_i-\overline{y})}}}$

(6)$RMSE=\sqrt{{\displaystyle \frac{{\displaystyle {\sum }_{i=1}^n{(y_i-{\widehat{y}}_i)}^2}}{n}}}$

(7)$rRMSE={\displaystyle \frac{RMSE}{\overline{y}}}\times 100\%$

where $y_i$ is the reference value of AGB, ${\widehat{y}}_i$ is the estimated value of AGB, $\overline{y}$ is the average of AGB measured values, and n is the number of samples.

Semivariance functions are evaluated using$R^2$and the residual sum of squares (RSS). The structural ratio (SR) is employed to determine the spatial correlation of the forest AGB. Subsequently, the optimal interpolation model is selected based on the criteria of having a larger$R^2$, a higher SR, and a smaller RSS. The formulas are as follows:

(8)$RSS={\displaystyle \sum _{i=1}^n{(s_i-{\widehat{s}}_i)}^2}$

where $s_i$ is the actual semivariance value, ${\widehat{s}}_i$ is the predicted semivariance value, and n is the number of values. A lower RSS value indicates relatively better estimation performance of the semivariance function model [68].

(9)$SR={\displaystyle \frac{C_1}{C_0+C_1}}$

where SR represents the structure ratio, $C_0$ is nugget, $C_1$ is partial sill, and $C_0+C_1$ is sill.

The accuracy of the co-kriging model is evaluated through cross-validation in the geostatistical analysis module [70]. Besides the RMSE, the mean error (ME), mean standardized error (MSE), root mean square standardized error (RMSSE), and average standard error (ASE) are also provided in this process to assess the performance of the interpolation model.

(10)$ME={\displaystyle \frac{{\displaystyle {\sum }_{i=1}^n(v_i-{\widehat{v}}_i)}}{n}}$

(11)$MSE={\displaystyle \frac{{\displaystyle {\sum }_{i=1}^n[({\widehat{v}}_i-v_i)/{\widehat{\sigma }}_i]}}{n}}$

(12)$RMSSE=\sqrt{{\displaystyle \frac{{\displaystyle {\sum }_{i=1}^n{[({\widehat{v}}_i-v_i)/{\widehat{\sigma }}_i]}^2}}{n}}}$

(13)$ASE=\sqrt{{\displaystyle \frac{{\displaystyle {\sum }_{i=1}^n{{\widehat{\sigma }}_i}^2}}{n}}}$

where $v_i$ represents the observed value, ${\widehat{v}}_i$ represents the predicted value, ${\widehat{\sigma }}_i$ is the standard error of predicted values, and n is the number of samples used.

4. Results and Analysis

4.1. Variable Selection Result

The Pearson correlation coefficient was used to analyze the correlation between AGB and all features. From Figure 4, we can see that the optimized features, including rh95, rh96, cover, fhd_normal, pai, B12, B4, MNDWI, NDVI, and B2_savg, have a strong correlation with biomass (0.802, 0.788, 0.621, 0.673, 0.588, −0.688, −0.644, 0.721, 0.649, and −0.633), all of which correlated with AGB at the 0.01 level. Among them, certain features exhibited a strong negative correlation with AGB. This may be attributed to the fact that as these feature values increase, the proportion of the area covered by buildings increases or vegetation growth deteriorates, leading to relatively lower AGB values [59]. The other features related to topography and vegetation were also correlated with the biomass in the study area, but the correlation is slightly weaker.

The stepwise method was employed to further optimize the independent variables of the MLR model. The results show that as the number of variables increases to 47, the AIC value tends to be minimized. Specifically, as fhd_normal, landsat_treecover, rv_a1, rh95, MNDWI, and B4_asm were added to the model gradually, the AIC values decreased rapidly (see Figure 5a). rh95 exhibited the most pronounced correlation, with the highest explanatory power for forest biomass, followed by landsat_treecover. For non-parametric feature selection, 45 variables were considered because their cumulative percentage of eigenvalues (CPE) reached 90% [71]. From Figure 5b, we can see that rh95 had the highest importance score, significantly higher than the other variables, followed by the B12 band and landsat_treecover. The results of the variable selection are similar to the stepwise method. Canopy cover and rh95 can reflect the structural characteristics and health status of the forest [72]. B12 is the short-wave infrared band, which is sensitive to minor changes in vegetation canopy structure and the growth status of vegetation [73]. Therefore, these variables have a significant impact on AGB prediction.

4.2. Footprint-Level AGB Estimation Results

In this study, the data were randomly divided into an 80% training set and a 20% validation set, considering the overall distribution of the data and ensuring the accuracy of the interpolated sample points (i.e., the footprint-level AGB). Among the three modeling schemes designed based on different data sources, the CatBoost model consistently outperformed the other models in terms of estimation precision. The RF model followed, while the predictive performance of the MLR model was the poorest. As can be seen from Table 3, the MLR model has the longest runtime and the most complex parameter settings (i.e., the highest number of variables). The performance of the RF model was close to that of CatBoost, but RF was inferior to CatBoost in terms of estimation accuracy and computational efficiency. The CatBoost model demonstrated relatively good performance thanks to the strategy of gradient boosting, which reduced the risk of overfitting by optimizing the gradient bias [74]. The MLR model is sensitive to outliers and has limitations on modeling in complex scenarios [42,75]. Therefore, the accuracy of this model was the lowest among the three modeling schemes.

From Table 3, we can see that when the single data source was used for the forest AGB inversion, the accuracy of GEDI data was consistently higher than that of Sentinel-2 data. For the three regression models, compared to Sentinel-2, the ${\mathrm{R}}^2$for GEDI predictions increased by 0.1, 0.1, and 0.08; the RMSE decreased by 15.87 Mg/ha, 15.63 Mg/ha, and 10.34 Mg/ha; and the rRMSE decreased by 8.54%, 8.62%, and 5.64%, respectively. When feature variables from GEDI and Sentinel-2 data were fused for the AGB inversion, the CatBoost model achieved the highest accuracy, with an ${\mathrm{R}}^2$ of 0.87, an RMSE of 49.56 Mg/ha, and an rRMSE of 27.06%. This result indicates that the combination of multi-source data can effectively mitigate the signal saturation associated with optical images, thus enhancing the accuracy and reliability of forest AGB estimation [76].

Figure 6 shows the regression results of the three models using features from GEDI and Sentinel-2 data. We found that all three models exhibited underestimation in the high-value range, which may be due to the uneven distribution of biomass. The number of samples in the high-value range is relatively smaller compared to the low-value range and includes some noise and outliers, thus affecting the models’ predictions for high AGB values. The outlier phenomenon in the MLR model was more significant compared to the other two models, with a bias of −0.75 Mg/ha; RF was second with a bias of −0.38 Mg/ha. The CatBoost model had the best fitting performance with a bias of −0.29 Mg/ha, showing higher consistency between the predicted values and the reference values. Therefore, CatBoost was employed in this study to establish a forest AGB prediction model at the footprint level using the fused feature variables. The prediction results of this model served as the main variable for co-kriging interpolation mapping.

4.3. Wall-to-Wall Forest AGB Mapping Results

4.3.1. Accuracy Analysis of Forest AGB Mapping Results

To verify the accuracy of co-kriging interpolation mapping, 77,371 GEDI footprints were randomly divided into a training set and a validation set at a ratio of 8:2. As a result, a total of 61,895 GEDI footprints were involved in the interpolation mapping, and the remaining randomly selected 15,476 footprints were used to assess the mapping precision. The results of interpolation mapping are shown in Table 4.

From Table 4, the cross-validation results of the interpolation model show that the MSE is near 0, the RMSSE is close to 1, and the value of ASE is close to the RMSE, which demonstrates that the predictions are unbiased. The mean AGB of the interpolated predictions in this study is close to that of the ALS data across a large number of samples, reflecting the unbiased nature of the interpolated results. However, the difference between the average (Ave) and the standard deviation (SD) Indicates a significant dispersion of AGB values in this area, reflecting an uneven distribution of AGB, with some areas exhibiting extreme values. From Figure 7, it can be seen that the AGB values are primarily distributed in the range of 0–50 Mg/ha, which is much higher than in the other ranges. Compared to the reference values from ALS, the predicted values from the co-kriging interpolation exhibit a broader distribution in the range of 50–200 Mg/ha. As biomass increases, the number of occurrences gradually decreases, with a few AGB values exceeding 700 Mg/ha. Overall, the AGB distribution derived from co-kriging aligns with the distribution of the ALS data, further demonstrating the unbiased nature of the interpolation predictions in this study.

In addition to performing cross-validation on the interpolation model, the accuracy of the wall-to-wall forest AGB map was also validated, with an ${\mathrm{R}}^2$ of 0.69, an RMSE of 81.56 Mg/ha, and an rRMSE of 40.98%. As shown in Figure 8, the distribution of AGB values shows that half of the AGB values fall between 100 and 300 Mg/ha, and the data overall exhibit a skewed distribution with a few extreme high values. The overall distribution of AGB derived from the interpolation method is almost consistent with that of the airborne AGB within the 1.5 IQR (interquartile range). Compared to the airborne reference values, the AGB obtained via co-kriging is relatively less distributed in high-value ranges, with a bias of −3.236 Mg/ha. The result indicates that the predictions tend to be underestimated. This may be due to the underestimation of the variability in sample points (i.e., where RMSE exceeds ASE and RMSSE is greater than 1). Another reason for this result may be the relatively small number of high-value samples, leading to the model’s inadequate response to extreme high values in biomass interpolation.

Overall, forest AGB values estimated in this study area primarily range between 0 and 600 Mg/ha. When the AGB value exceeds 600 Mg/ha, the interpolation method can still be used to fit biomass, but the predicted values are underestimated to some extent. When the AGB value is less than 600 Mg/ha within the normal range, the accuracy of the spatially continuous forest AGB mapped using co-kriging is relatively high, with lower error. These results indicate that the interpolation method has potential for mapping high-quality spatially continuous forest AGB.

4.3.2. Spatial Distribution Characteristics of Wall-to-Wall Forest AGB

Figure 9 displays the remote-sensing images of the study area and the mapping results of the forest AGB at a 30 m resolution. Non-vegetated areas within the AGB product were removed using the Global 30 m Fine Ground Cover 2020 product. Forest AGB values are higher in the north-western part of Sonoma County compared to the eastern and southern areas. The north-western part of the study area, which is bordered by the ocean, receives ample precipitation that supports healthy vegetation. Furthermore, significant topographical variations and minimal human activity contribute to generally higher AGB values in this area. In contrast, lower forest AGB values are observed in the eastern and southern areas, which are close to the central residential parts of the county. These residential areas are characterized by high levels of human activity, and the existing vegetation has been destroyed by the construction of towns and roads, resulting in lower AGB values. Artifacts were observed in areas with fewer GEDI footprints due to the uneven distribution of these footprints, but the strip effect inherent in the GEDI data itself was mitigated. As shown in Figure 9, the spatial distribution of the forest AGB is generally consistent with the vegetation distribution in the study area.

5. Discussion

5.1. Selection of the Interpolation Model

Before interpolation, it is necessary to determine the optimal semivariogram model that fits the correlation among the AGB sample points. The precision of the forest AGB mapping can be improved by considering this correlation during the interpolation process [77,78]. In this study, the Gauss, spherical, linear, and exponential functions were compared to determine the optimal model for footprint-level AGB interpolation. The fitting results of each model are shown in Table 5.

From Table 5, we can see that the exponential model achieved relatively higher ${\mathrm{R}}^2$ values and lower RSS values for the main variable (i.e., footprint-level predicted AGB) and the covariates (i.e., rh95, B12, and slope). The SR is used to measure the degree of spatial autocorrelation of the variables [36]. If the SR > 0.75, the variable has strong spatial correlation; if the SR < 0.25, the spatial correlation is weak; if the SR lies between 0.25 and 0.75, the spatial correlation is moderately strong [79]. The SR of the predicted AGB reached 0.82 when using the exponential model to fit the correlation among the AGB sample points, indicating a strong spatial correlation of footprint-level AGB within the study area. The result that the exponential model is more applicable in this study reveals that the distribution of AGB in the area exhibits a short-distance correlation. Specifically, the spatial correlation of AGB rapidly weakens beyond a certain distance within the region [80], which may be due to significant local environmental differences, such as soil types, moisture conditions, and microclimatic changes. Consequently, compared with the default model for interpolation, the interpolation results obtained from the optimal model selected by calculating the semivariance function have a higher accuracy [42].

5.2. The Impact of Covariate Selection on the Interpolation Accuracy of AGB

Co-kriging extends the optimal estimation of regionalized variables from a single attribute to two or more collaborating regionalized attributes [81]. It is crucial to select the appropriate number of covariates that are highly correlated with the main variable, because this step can effectively improve the performance and accuracy of co-kriging interpolation [82]. In this study, four experiments were designed to compare the effects of different covariate combinations on interpolation results. The combinations are as follows: no covariates; only the rh95 from GEDI; rh95 and Sentinel-2 B12; and rh95, B12, and slope, which can be seen in Table 6. In the experiments for interpolation mapping, the footprint-level predicted AGB served as the main variable, while rh95, B12, and slope were selected as covariates based on the integrated results of correlation analysis, feature importance analysis, and the influence of terrain.

Table 6 shows the interpolation accuracy using different covariates. The interpolation accuracy of ordinary kriging (no covariates) was the lowest, with an ${\mathrm{R}}^2$ value of 0.62, an RMSE of 91.06 Mg/ha, and an rRMSE of 45.76%. In comparison, the inclusion of covariates effectively improves interpolation accuracy and reduces prediction errors, as the interpolation process fully considers the statistical correlations and spatial relationships between the variables [83]. When the number of covariates was two—specifically, rh95 and B12 were selected as covariates—the interpolation accuracy was the highest, with an ${\mathrm{R}}^2$value of 0.69, an RMSE of 81.56 Mg/ha, and an rRMSE of 40.98%. The results indicate that tree height and B12 can provide additional useful information about forest status and AGB distribution, thus improving the accuracy of spatial interpolation. By increasing the number of covariates, it is possible to provide more information, which theoretically helps in more accurately estimating the spatial distribution of the main variable. However, if too many covariates are involved in co-kriging interpolation, it may lead to overfitting of the model and reduce the generalizability, especially in the case of uneven distribution of sampling points [84]. From Figure 10, we can see that the fitting accuracy begins to decline and the prediction error gradually increases when the number of covariates reaches two. As a result, the rh95 and B12 were selected as the covariates for the interpolation of forest AGB in this study, because the best result was achieved by using this combination.

5.3. Accuracy Analysis of Interpolation and Regression Mapping Results

In traditional methods, optical images are commonly employed to map the spatial distribution of forest biomass [29,65,85]. Specifically, the features selected by parametric or non-parametric methods are used as inputs in machine-learning algorithms to build biomass extrapolation models. However, machine-learning regression methods based solely on optical remote-sensing images are susceptible to signal saturation problems, leading to significant prediction bias [86]. To mitigate the saturation effect from optical images, the co-kriging interpolation method was employed in this study to map forest AGB across the study area. We used 15,476 randomly selected footprints as an independent validation set and employed airborne-derived AGB values as reference AGB to evaluate the mapping accuracy. The accuracy comparison between this method and the regression method based on optical images (i.e., Sentinel-2A in this study) is shown in Figure 11.

From Figure 11, it can be seen that the mapping results of interpolation are better than those of the regression method. The ${\mathrm{R}}^2$value increased from 0.60 to 0.69, the RMSE decreased from 90.07 to 81.56 Mg/ha, and the rRMSE decreased by 3.83%. In areas with higher biomass, AGB values based on optical image regression are severely underestimated, while the performance of the interpolation method is relatively better. The bias of the regression results is −7.309 Mg/ha, while the bias of the interpolation results is −3.236 Mg/ha. Compared to the regression method, the prediction bias of the interpolation method decreased by 55%, with a significant decrease observed particularly when the AGB value exceeds 600 Mg/ha. These results demonstrate that the strategy of regression inversion does not take full advantage of the structural information from spaceborne LiDAR data. The accurate forest AGB information at GEDI footprints may be blurred during the regression process [43]. For interpolation techniques, the information from unsaturated areas can be utilized to influence and adjust estimation results in saturated areas by considering the relationships between geographical locations, thus minimizing errors caused by local saturation [81]. As a result, the AGB product generated by extrapolation based on the optical images is still affected by saturation effects in areas with higher forest AGB. Interpolated mapping can mitigate these saturation effects while improving mapping accuracy.

6. Conclusions

In this study, we proposed a co-kriging-interpolation-based method to map spatially continuous forest AGB by integrating GEDI and Sentinel-2 data. Specifically, we draw the following conclusions: (1) For the inversion of the footprint-level AGB, the CatBoost model produced the best estimation results compared to the other machine-learning algorithms (e.g., RF and MLR). (2) Compared to using optical images or LiDAR data alone, the combination of GEDI and Sentinel-2 data provided more accurate predictions for estimating footprint-level AGB. (3) Compared with the optical imagery extrapolation-based (i.e., regression) forest AGB mapping approach, co-kriging interpolation mitigated the saturation effect in areas with higher forest AGB and improved mapping accuracy, with a decrease in bias to −3.236 Mg/ha and an increase in ${\mathrm{R}}^2$to 0.69. Overall, the combination of spaceborne LiDAR and multispectral images enhances the estimation accuracy of forest AGB at the footprint level. Additionally, mapping accuracy is further improved by using the co-kriging method, which demonstrates great potential for monitoring large-scale forest biomass.

There are still some limitations to this study. For example, artifacts appeared in the interpolation results in areas with sparse sample data. This occurred because the interpolation algorithm may produce discontinuous or unnatural transitions when attempting to fill in large areas of unknown data. In future research, we aim to address this issue by integrating other spaceborne LiDAR data (e.g., ICESat-2/ATLAS) to increase the density of sampling points and improve the spatially continuous mapping of forest AGB.

Author Contributions

Conceptualization, Y.W. and H.W.; methodology, Y.W. and H.W.; software, C.W. and S.Z.; validation, Y.W.; investigation, R.W. and S.W.; resources, J.D.; data curation, Y.W.; writing—original draft preparation, Y.W.; writing—review and editing, Y.W. and H.W.; visualization, Y.W.; supervision, C.W.; project administration, H.W.; funding acquisition, H.W. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

GEDI data used in this study are available at https://search.earthdata.nasa.gov/ (accessed on 15 August 2023); airborne LiDAR data used in this study are available at https://daac.ornl.gov/ (accessed on 10 August 2023); Sentinel-2 data used in this study are openly available on the Google Earth Engine platform.

Acknowledgments

We would like to express our special thanks to NASA for providing the GEDI and airborne data, and to the editors and reviewers for their valuable comments.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Location of the study area. (a) The administrative boundary of the State of California and the location of the study area within the state (i.e., where the red box is); (b) ground elevation distribution of the study area; (c) local distribution of GEDI footprints.

Figure 2. Technology road map of this study.

Figure 3. Filtered distribution of GEDI footprints: (a) urban area; (b) vegetation area. The green area represents vegetation, the purple area represents urban areas, and the blue area represents the ocean.

Figure 4. The correlation between the features and forest AGB, with all significance levels of the selected features at 0.01.

Figure 5. Feature selection outcomes. (a) Features selected using the stepwise method. (b) Features selected using the random forest method.

Figure 6. Scatter plots between reference AGB and predicted AGB derived from the combined optical and GEDI data. The red dotted lines represent the fitted trend-lines.

Figure 7. Histograms of AGB: (a) ALS-derived and (b) interpolation-method-derived.

Figure 8. Boxplot of AGB distribution. IQR = Q3 − Q1, where Q1 is the first quartile and Q3 is the third quartile.

Figure 9. Forest AGB mapping results. (a) Satellite images of the study area. (b) Forest AGB map derived from the co-kriging interpolation model. The red line represents the boundary of the study area (i.e., Sonoma County).

Figure 10. Comparison of interpolation accuracy for different covariate combinations.

Figure 11. Accuracy assessment of the wall-to-wall forest AGB products. (a) The predicted AGB derived from co-kriging interpolation. (b) The predicted AGB derived from the regression method. The red dotted lines represent the fitted trend-lines.

Table 1

Description of the GEDI L2 products used in this study.

Product	Parameters	Level	Resolution
L2A	Canopy height metrics	Footprint level	25 m
L2B	Canopy cover(total and vertical profiles of canopy cover),	Footprint level	25 m
	Plant Area Index (PAI),
	Plant Area Volume Density (PAVD),
	Foliage Height Diversity (FHD)

Table 2

Features extracted from GEDI, Sentinel-2, and SRTM DEM.

Source	Type/Number	Feature	Description
GEDI	L2A (11)	rh90, rh91, rh92, rh93, rh94, rh95, rh96, rh97, rh98, rh99, rh100	Relative height metrics
	L2B (28)	cover, pai, fhd_normal, paga_theta, landsat_treecover,	The footprint-level canopy coverage and vertical profile metrics
		leaf_off_doy, leaf_on_doy, leaf_off_flag, modis_treecover,
		modis_nonvegetated, rg_aN, rv_aN, rx_energy_aN
Sentinel-2	Spectral band (10)	B2, B3, B4, B5, B6, B7, B8, B8a, B11, B12	Band reflectance information
	Spectral index (8)	NDVI, normalized difference vegetation index	(B8 − B4)/(B8 + B4)
		RVI, relative vegetation index	B8/B4
		EVI, enhanced vegetation index,	2.5 × [(B8 − B4)/ (B8 + 6 × B4 − 7.5 × B2 + 1)]
		DVI, difference vegetation index	B8 − B4
		MNDWI, modified normalized difference water index	(B3 − B11)/(B3 + B11)
		RGVI, red–green vegetation index	(B4 − B3)/(B4 + B3)
		SAVI, soil-adjusted vegetation index	[(B8 − B4)/(B8 + B4 + L)] × (1 + L)
		SVI, shadow vegetation index	(B8 − B4) × B8/(B8 + B4)
	Textural feature (170)	ASMi, CONTRASTi, CORRi, VARi, IDMi, SAVGi, SVAGi,	Reflecting surface characteristics of land cover
		SENTi, ENTi, DVARi, DENTi, IMCORR1i, IMCORR2i,
		DISSi, INERTIAi, SHADEi, PROMi
SRTM DEM	Terrain factor (3)	Elevation	Topographic-feature-related factors
		Slope
		Aspect

L is 0.5; _aN (N = 1~6), which means 6 algorithms for GEDI; i represents the 10 bands (B2, B3, B4, B5, B6, B7, B8, B8a, B11, and B12) extracted from the Sentinel-2 image.

Table 3

Performance comparison of different models and data sources.

Data	Model	Features	Runtime (s)	${\mathbf{R}}^2$	RMSE (Mg/ha)	rRMSE (%)
GEDI	CatBoost	16	45	0.84	55.32	30.23
	RF	16	62	0.82	59.49	32.07
	MLR	31	85	0.77	67.71	36.69
Sentinel-2	CatBoost	65	82	0.74	71.19	38.77
	RF	65	98	0.72	75.12	40.69
	MLR	94	145	0.69	78.05	42.33
GEDI + Sentinel-2	CatBoost	45	65	0.87	49.56	27.06
	RF	45	78	0.85	54.08	29.37
	MLR	47	117	0.81	60.69	32.91

Table 4

Accuracy of co-kriging interpolation with the exponential model.

Cross-Validation Accuracy	ME (Mg/ha)	RMSE (Mg/ha)	MSE (Mg/ha)2	ASE (Mg/ha)2	RMSSE
	0.99	65.20	0.01	61.82	1.09
Mapping Accuracy	${\mathrm{R}}^2$	RMSE (Mg/ha)	rRMSE (%)	BIAS (Mg/ha)
	0.69	81.56	40.98	−3.236

Table 5

Correlation fitting results with different models.

	Variable	Model	Nugget	Sill	SR	Range (m)	${\mathbf{R}}^2$	RSS
Main variable	Predicted AGB	Linear	0.013	0.019	0.27	56,837.21	0.73	7.44 × 10−6
		Spherical	0.001	0.016	0.97	6800.00	0.49	1.40 × 10−6
		Exponential	0.008	0.045	0.82	31,800.00	0.93	2.55 × 10−6
		Gaussian	0.002	0.016	0.87	5715.76	0.49	1.40 × 10−6
Covariates	rh95	Linear	1.445	1.980	0.27	57,383.05	0.94	0.0181
		Spherical	0.053	1.753	0.97	5200.00	0.33	0.215
		Exponential	1.408	2.817	0.50	339,300.00	0.97	0.0103
		Gaussian	0.221	1.753	0.87	4503.33	0.33	0.215
	B12	Linear	0.276	0.402	0.31	57,383.05	0.90	1.50 × 10−4
		Spherical	0.006	0.352	0.98	6800.00	0.45	8.67 × 10−4
		Exponential	0.017	0.147	0.88	237,300.00	0.93	9.79 × 10−4
		Gaussian	0.039	0.352	0.89	5888.97	0.45	8.65 × 10−4
	Slope	Linear	0.011	0.014	0.21	57,383.05	0.86	1.62 × 10−6
		Spherical	0.001	0.013	0.96	5100.00	0.44	6.60 × 10−6
		Exponential	0.002	0.013	0.87	6900.00	0.54	5.57 × 10−6
		Gaussian	0.001	0.013	0.87	4330.12	0.44	6.59 × 10−6

Table 6

Interpolation results of different covariate combinations.

Covariate	${\mathbf{R}}^2$	RMSE (Mg/ha)	rRMSE (%)	BIAS (Mg/ha)
-	0.62	91.06	45.76	−4.474
rh95	0.65	87.35	43.89	−3.024
rh95 + B12	0.69	81.56	40.98	−3.236
rh95 + B12 + slope	0.66	85.83	43.13	−3.231

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/rs16162913/s1, Table S1: Features extracted from GEDI; Table S2: Features extracted from Sentinel-2; Table S3: 170 texture features; Table S4: Features extracted from SRTM DEM.

References

1. FAO. The State of the World’s Forests 2018: Forest Pathways to Sustainable Development; UN: New York, NY, USA, 2018.

2. Zhang, C.; Ju, W.; Chen, J.M.; Wang, X.; Yang, L.; Zheng, G. Disturbance-induced reduction of biomass carbon sinks of China’s forests in recent years. Environ. Res. Lett.; 2015; 10, 114021. [DOI: https://dx.doi.org/10.1088/1748-9326/10/11/114021]

3. Le Toan, T.; Quegan, S.; Davidson, M.; Balzter, H.; Paillou, P.; Papathanassiou, K.; Plummer, S.; Rocca, F.; Saatchi, S.; Shugart, H. The BIOMASS mission: Mapping global forest biomass to better understand the terrestrial carbon cycle. Remote Sens. Environ.; 2011; 115, pp. 2850-2860. [DOI: https://dx.doi.org/10.1016/j.rse.2011.03.020]

4. Wang, P.; Tan, S.; Zhang, G.; Wang, S.; Wu, X. Remote sensing estimation of forest aboveground biomass based on Lasso-SVR. Forests; 2022; 13, 1597. [DOI: https://dx.doi.org/10.3390/f13101597]

5. Chen, Q.; McRoberts, R.E.; Wang, C.; Radtke, P.J. Forest aboveground biomass mapping and estimation across multiple spatial scales using model-based inference. Remote Sens. Environ.; 2016; 184, pp. 350-360. [DOI: https://dx.doi.org/10.1016/j.rse.2016.07.023]

6. Sha, Z.; Bai, Y.; Li, R.; Lan, H.; Zhang, X.; Li, J.; Liu, X.; Chang, S.; Xie, Y. The global carbon sink potential of terrestrial vegetation can be increased substantially by optimal land management. Commun. Earth Environ.; 2022; 3, 8. [DOI: https://dx.doi.org/10.1038/s43247-021-00333-1]

7. Liang, X.; Hyyppä, J.; Kaartinen, H.; Lehtomäki, M.; Pyörälä, J.; Pfeifer, N.; Holopainen, M.; Brolly, G.; Francesco, P.; Hackenberg, J. et al. International benchmarking of terrestrial laser scanning approaches for forest inventories. ISPRS J. Photogramm. Remote Sens.; 2018; 144, pp. 137-179. [DOI: https://dx.doi.org/10.1016/j.isprsjprs.2018.06.021]

8. Wang, Y.; Lehtomäki, M.; Liang, X.; Pyörälä, J.; Kukko, A.; Jaakkola, A.; Liu, J.; Feng, Z.; Chen, R.; Hyyppä, J. et al. Is field-measured tree height as reliable as believed–A comparison study of tree height estimates from field measurement, airborne laser scanning and terrestrial laser scanning in a boreal forest. ISPRS J. Photogramm. Remote Sens.; 2019; 147, pp. 132-145. [DOI: https://dx.doi.org/10.1016/j.isprsjprs.2018.11.008]

9. Arévalo, P.; Baccini, A.; Woodcock, C.E.; Olofsson, P.; Walker, W.S. Continuous mapping of aboveground biomass using Landsat time series. Remote Sens. Environ.; 2023; 288, 113483. [DOI: https://dx.doi.org/10.1016/j.rse.2023.113483]

10. David, R.M.; Rosser, N.J.; Donoghue, D.N. Improving above ground biomass estimates of Southern Africa dryland forests by combining Sentinel-1 SAR and Sentinel-2 multispectral imagery. Remote Sens. Environ.; 2022; 282, 113232. [DOI: https://dx.doi.org/10.1016/j.rse.2022.113232]

11. Qin, S.; Nie, S.; Guan, Y.; Zhang, D.; Wang, C.; Zhang, X. Forest emissions reduction assessment using airborne LiDAR for biomass estimation. Resour. Conserv. Recycl.; 2022; 181, 106224. [DOI: https://dx.doi.org/10.1016/j.resconrec.2022.106224]

12. Silveira, E.M.; Radeloff, V.C.; Martinuzzi, S.; Pastur, G.J.M.; Bono, J.; Politi, N.; Lizarraga, L.; Rivera, L.O.; Ciuffoli, L.; Rosas, Y.M. Nationwide native forest structure maps for Argentina based on forest inventory data, SAR Sentinel-1 and vegetation metrics from Sentinel-2 imagery. Remote Sens. Environ.; 2023; 285, 113391. [DOI: https://dx.doi.org/10.1016/j.rse.2022.113391]

13. Lechner, A.M.; Foody, G.M.; Boyd, D.S. Applications in remote sensing to forest ecology and management. One Earth; 2020; 2, pp. 405-412. [DOI: https://dx.doi.org/10.1016/j.oneear.2020.05.001]

14. Moradi, F.; Sadeghi, S.M.M.; Heidarlou, H.B.; Deljouei, A.; Boshkar, E.; Borz, S.A. Above-ground biomass estimation in a Mediterranean sparse coppice oak forest using Sentinel-2 data. Ann. For. Res.; 2022; 65, pp. 165-182. [DOI: https://dx.doi.org/10.15287/afr.2022.2390]

15. Qian, C.; Qiang, H.; Wang, F.; Li, M. Estimation of Forest Aboveground Biomass in Karst Areas Using Multi-Source Remote Sensing Data and the K-DBN Algorithm. Remote Sens.; 2021; 13, 5030. [DOI: https://dx.doi.org/10.3390/rs13245030]

16. López-Serrano, P.M.; López Sánchez, C.A.; Solís-Moreno, R.; Corral-Rivas, J.J. Geospatial estimation of above ground forest biomass in the Sierra Madre Occidental in the state of Durango, Mexico. Forests; 2016; 7, 70. [DOI: https://dx.doi.org/10.3390/f7030070]

17. Kumar, P.; Krishna, A.P. Forest biomass estimation using multi-polarization SAR data coupled with optical data. Curr. Sci.; 2020; 119, pp. 1316-1321. [DOI: https://dx.doi.org/10.18520/cs/v119/i8/1316-1321]

18. Zhang, W.; Zhang, Y.; Fan, W.; Wu, G. Comparison of the accuracy of forest biomass estimation by interference water cloud model for sentinel data with different polarization modes. J. Northeast For. Univ.; 2020; 48, pp. 27-32.

19. Ji, Y.; Wang, L.; Zhang, W.; Marino, A.; Wang, M.; Ma, J.; Shi, J.; Jing, Q.; Zhang, F.; Zhao, H. Forest above-ground biomass estimation using X, C, L, and P band SAR polarimetric observations and different inversion models. Int. J. Digit. Earth; 2024; 17, 2310730. [DOI: https://dx.doi.org/10.1080/17538947.2024.2310730]

20. Guo, Q.; Su, Y.; Hu, T.; Guan, H.; Jin, S.; Zhang, J.; Zhao, X.; Xu, K.; Wei, D.; Kelly, M. et al. Lidar boosts 3D ecological observations and modelings: A review and perspective. IEEE Geosci. Remote Sens. Mag.; 2020; 9, pp. 232-257. [DOI: https://dx.doi.org/10.1109/MGRS.2020.3032713]

21. Beland, M.; Parker, G.; Sparrow, B.; Harding, D.; Chasmer, L.; Phinn, S.; Antonarakis, A.; Strahler, A. Management. On promoting the use of lidar systems in forest ecosystem research. For. Ecol. Manag.; 2019; 450, 117484. [DOI: https://dx.doi.org/10.1016/j.foreco.2019.117484]

22. Kozak, I.; Popov, M.; Semko, I.; Mylenka, M.; Kozak-Balaniuk, I.J.U.F.; Greening, U. Improving methods to predict aboveground biomass of Pinus sylvestris in urban forest using UFB model, LiDAR and digital hemispherical photography. Urban For. Urban Green.; 2023; 79, 127793. [DOI: https://dx.doi.org/10.1016/j.ufug.2022.127793]

23. Zhang, L.; Zhao, Y.; Chen, C.; Li, X.; Mao, F.; Lv, L.; Yu, J.; Song, M.; Huang, L.; Chen, J. UAV-LiDAR Integration with Sentinel-2 Enhances Precision in AGB Estimation for Bamboo Forests. Remote Sens.; 2024; 16, 705. [DOI: https://dx.doi.org/10.3390/rs16040705]

24. Hartley, R.J.; Leonardo, E.M.; Massam, P.; Watt, M.S.; Estarija, H.J.; Wright, L.; Melia, N.; Pearse, G.D. An assessment of high-density UAV point clouds for the measurement of young forestry trials. Remote Sens.; 2020; 12, 4039. [DOI: https://dx.doi.org/10.3390/rs12244039]

25. Padalia, H.; Prakash, A.; Watham, T.J. Modelling aboveground biomass of a multistage managed forest through synergistic use of Landsat-OLI, ALOS-2 L-band SAR and GEDI metrics. Ecol. Inf.; 2023; 77, 102234. [DOI: https://dx.doi.org/10.1016/j.ecoinf.2023.102234]

26. Pang, Y.; Li, Z.; Chen, B.; Liang, X. Progress and trend of spaceborne lidar forest detection. Aerosp. Shanghai; 2019; 36, pp. 20-28.

27. Salas, E.A.L. Waveform LiDAR concepts and applications for potential vegetation phenology monitoring and modeling: A comprehensive review. Geo-Spat. Inf. Sci.; 2021; 24, pp. 179-200. [DOI: https://dx.doi.org/10.1080/10095020.2020.1761763]

28. Varvia, P.; Saarela, S.; Maltamo, M.; Packalen, P.; Gobakken, T.; Næsset, E.; Ståhl, G.; Korhonen, L. Estimation of boreal forest biomass from ICESat-2 data using hierarchical hybrid inference. Remote Sens. Environ.; 2024; 311, 114249. [DOI: https://dx.doi.org/10.1016/j.rse.2024.114249]

29. Liang, M.; Duncanson, L.; Silva, J.A.; Sedano, F. Quantifying aboveground biomass dynamics from charcoal degradation in Mozambique using GEDI Lidar and Landsat. Remote Sens. Environ.; 2023; 284, 113367. [DOI: https://dx.doi.org/10.1016/j.rse.2022.113367]

30. Shendryk, Y. Fusing GEDI with earth observation data for large area aboveground biomass mapping. Int. J. Appl. Earth Obs. Geoinf.; 2022; 115, 103108. [DOI: https://dx.doi.org/10.1016/j.jag.2022.103108]

31. Chen, L.; Ren, C.; Bao, G.; Zhang, B.; Wang, Z.; Liu, M.; Man, W.; Liu, J. Improved object-based estimation of forest aboveground biomass by integrating LiDAR data from GEDI and ICESat-2 with multi-sensor images in a heterogeneous mountainous region. Remote Sens.; 2022; 14, 2743. [DOI: https://dx.doi.org/10.3390/rs14122743]

32. Narine, L.L.; Popescu, S.C.; Malambo, L. Using ICESat-2 to estimate and map forest aboveground biomass: A first example. Remote Sens.; 2020; 12, 1824. [DOI: https://dx.doi.org/10.3390/rs12111824]

33. Wang, C.; Zhang, W.; Ji, Y.; Marino, A.; Li, C.; Wang, L.; Zhao, H.; Wang, M. Estimation of Aboveground Biomass for Different Forest Types Using Data from Sentinel-1, Sentinel-2, ALOS PALSAR-2, and GEDI. Forests; 2024; 15, 215. [DOI: https://dx.doi.org/10.3390/f15010215]

34. Wang, Y.; Peng, Y.; Hu, X.; Zhang, P. Fine-Resolution Forest Height Estimation by Integrating ICESat-2 and Landsat 8 OLI Data with a Spatial Downscaling Method for Aboveground Biomass Quantification. Forests; 2023; 14, 1414. [DOI: https://dx.doi.org/10.3390/f14071414]

35. Nandy, S.; Srinet, R.; Padalia, H. Mapping forest height and aboveground biomass by integrating ICESat-2, Sentinel-1 and Sentinel-2 data using Random Forest algorithm in northwest Himalayan foothills of India. Geophys. Res. Lett.; 2021; 48, e2021GL093799. [DOI: https://dx.doi.org/10.1029/2021GL093799]

36. Song, H.; Xi, L.; Shu, Q.; Wei, Z.; Qiu, S. Estimate forest aboveground biomass of mountain by ICESat-2/ATLAS data interacting cokriging. Forests; 2022; 14, 13. [DOI: https://dx.doi.org/10.3390/f14010013]

37. Feng, Y. Spatial Statistics Theory and Its Application in Forestry; China Forestry Publishing House: Beijing, China, 2008.

38. Chiles, J.-P.; Delfiner, P. Geostatistics: Modeling Spatial Uncertainty; John Wiley & Sons: Hoboken, NJ, USA, 2012; Volume 713.

39. Zhao, M.; Zhao, Y.; Shen, T.; Li, S.; Yao, G.; Chen, Z.; Liu, Y.; Xu, Y. Research on the distribution of soil manganese and zinc pollution based on in-situ PXRF data, combined with Kriging interpolation and high-resolution mapping. Res. Environ. Sci.; 2023; 36, pp. 599-609.

40. Wang, C.; Qu, Y.; Shuai, Y. The collaborative Kriging method based on ArcGIS is used to predict the spatial distribution of gas pipeline network. Geomat. Spat. Inf. Technol.; 2022; 45, pp. 28-32.

41. Wang, H.; Zhao, M.; Huang, X.; Song, X.; Cai, B.; Tang, R.; Sun, J.; Han, Z.; Yang, J.; Liu, Y. Improving prediction of soil heavy metal (loid) concentration by developing a combined Co-kriging and geographically and temporally weighted regression (GTWR) model. J. Hazard. Mater.; 2024; 468, 133745. [DOI: https://dx.doi.org/10.1016/j.jhazmat.2024.133745]

42. Xu, L.; Shu, Q.; Fu, H.; Zhou, W.; Luo, S.; Gao, Y.; Yu, J.; Guo, C.; Yang, Z.; Xiao, J. Estimation of Quercus biomass in Shangri-La based on GEDI spaceborne LiDAR data. Forests; 2023; 14, 876. [DOI: https://dx.doi.org/10.3390/f14050876]

43. Liu, X.; Su, Y.; Hu, T.; Yang, Q.; Liu, B.; Deng, Y.; Tang, H.; Tang, Z.; Fang, J.; Guo, Q. Neural network guided interpolation for mapping canopy height of China’s forests by integrating GEDI and ICESat-2 data. Remote Sens. Environ.; 2022; 269, 112844. [DOI: https://dx.doi.org/10.1016/j.rse.2021.112844]

44. Li, H.; Nie, S.; Xi, X.; Wang, H.; Zhang, H.; Wang, C. Forest aboveground biomass inversion based on ICESat-2 / ATLAS and characteristic parameters. Laser J.; 2023; 44, pp. 62-67.

45. Silva, C.A.; Duncanson, L.; Hancock, S.; Neuenschwander, A.; Thomas, N.; Hofton, M.; Fatoyinbo, L.; Simard, M.; Marshak, C.Z.; Armston, J. Fusing simulated GEDI, ICESat-2 and NISAR data for regional aboveground biomass mapping. Remote Sens. Environ.; 2021; 253, 112234. [DOI: https://dx.doi.org/10.1016/j.rse.2020.112234]

46. Cooper, S.; Okujeni, A.; Pflugmacher, D.; van der Linden, S.; Hostert, P. Combining simulated hyperspectral EnMAP and Landsat time series for forest aboveground biomass mapping. Int. J. Appl. Earth Obs. Geoinf.; 2021; 98, 102307. [DOI: https://dx.doi.org/10.1016/j.jag.2021.102307]

47. CMS: LiDAR-Derived Biomass, Canopy Height and Cover, Sonoma County, California. 2013; Available online: https://daac.ornl.gov/CMS/guides/CMS_LiDAR_Biomass_CanHt_Sonoma.html (accessed on 10 August 2023).

48. Dubayah, R.; Blair, J.B.; Goetz, S.; Fatoyinbo, L.; Hansen, M.; Healey, S.; Hofton, M.; Hurtt, G.; Kellner, J.; Luthcke, S. The Global Ecosystem Dynamics Investigation: High-resolution laser ranging of the Earth’s forests and topography. Sci. Remote Sens.; 2020; 1, 100002. [DOI: https://dx.doi.org/10.1016/j.srs.2020.100002]

49. Wang, C. Accuracy Analysis and Improvement Methods for Forest Structure and Functioning Parameters of GEDI Products. Ph.D. Thesis; China University of Mining and Technology: Xuzhou, China, 2023.

50. Duan, J.; Wang, H.; Wang, C.; Nie, S.; Yang, X.; Xi, X. Denoising and classification of urban ICESat-2 photon data fused with Sentinel-2 spectral images. Int. J. Digit. Earth; 2023; 16, pp. 4346-4367. [DOI: https://dx.doi.org/10.1080/17538947.2023.2270513]

51. Zhou, Z. Research on Mangrove Remote Sensing Information Recognition Based on Multi-Source Remote Sensing Data. Master’s Thesis; Jilin University: Jilin, China, 2019.

52. Farr, T.G.; Rosen, P.A.; Caro, E.; Crippen, R.; Duren, R.; Hensley, S.; Kobrick, M.; Paller, M.; Rodriguez, E.; Roth, L. The shuttle radar topography mission. Rev. Geophys.; 2007; 45, [DOI: https://dx.doi.org/10.1029/2005RG000183]

53. Liu, Y.; Gong, W.; Xing, Y.; Hu, X.; Gong, J. Estimation of the forest stand mean height and aboveground biomass in Northeast China using SAR Sentinel-1B, multispectral Sentinel-2A, and DEM imagery. ISPRS J. Photogramm. Remote Sens.; 2019; 151, pp. 277-289. [DOI: https://dx.doi.org/10.1016/j.isprsjprs.2019.03.016]

54. 2020 Global 30 m Surface Coverage Fine Classification Product V1.0. Available online: https://data.casearth.cn/sdo/detail/5fbc7904819aec1ea2dd7061 (accessed on 12 February 2024).

55. Ngo, Y.-N.; Ho Tong Minh, D.; Baghdadi, N.; Fayad, I. Tropical forest top height by GEDI: From sparse coverage to continuous data. Remote Sens.; 2023; 15, 975. [DOI: https://dx.doi.org/10.3390/rs15040975]

56. Potapov, P.; Li, X.; Hernandez-Serna, A.; Tyukavina, A.; Hansen, M.C.; Kommareddy, A.; Pickens, A.; Turubanova, S.; Tang, H.; Silva, C.E. Mapping global forest canopy height through integration of GEDI and Landsat data. Remote Sens. Environ.; 2021; 253, 112165. [DOI: https://dx.doi.org/10.1016/j.rse.2020.112165]

57. Zeng, Y.; Hao, D.; Huete, A.; Dechant, B.; Berry, J.; Chen, J.M.; Joiner, J.; Frankenberg, C.; Bond-Lamberty, B.; Ryu, Y. et al. Optical vegetation indices for monitoring terrestrial ecosystems globally. Nat. Rev. Earth Environ.; 2022; 3, pp. 477-493. [DOI: https://dx.doi.org/10.1038/s43017-022-00298-5]

58. Yan, X.; Li, J.; Smith, A.R.; Yang, D.; Ma, T.; Su, Y.; Shao, J. Evaluation of machine learning methods and multi-source remote sensing data combinations to construct forest above-ground biomass models. Int. J. Digit. Earth; 2023; 16, pp. 4471-4491. [DOI: https://dx.doi.org/10.1080/17538947.2023.2270459]

59. Guo, Q.; Du, S.; Jiang, J.; Guo, W.; Zhao, H.; Yan, X.; Zhao, Y.; Xiao, W. Combining GEDI and sentinel data to estimate forest canopy mean height and aboveground biomass. Ecol. Inf.; 2023; 78, 102348. [DOI: https://dx.doi.org/10.1016/j.ecoinf.2023.102348]

60. Zhang, K.; Wu, X.; Zhao, J.; Liu, L.; Qin, P.; Wang, H.; Zhan, T. Real-time mountain fire risk assessment model of transmission corridor based on feature engineering, ensemble learning and model fusion. Power Syst. Technol.; 2023; 47, pp. 4727-4738.

61. Saber, M.; Boulmaiz, T.; Guermoui, M.; Abdrabo, K.I.; Kantoush, S.A.; Sumi, T.; Boutaghane, H.; Hori, T.; Binh, D.V.; Nguyen, B.Q. et al. Enhancing flood risk assessment through integration of ensemble learning approaches and physical-based hydrological modeling. Geomat. Nat. Hazards Risk; 2023; 14, 2203798. [DOI: https://dx.doi.org/10.1080/19475705.2023.2203798]

62. Xu, L.; Lai, H.; Yu, J.; Luo, S.; Guo, C.; Gao, Y.; Zhou, W.; Wang, S.; Shu, Q. Carbon Storage Estimation of Quercus aquifolioides Based on GEDI Spaceborne LiDAR Data and Landsat 9 Images in Shangri-La. Sustainability; 2023; 15, 11525. [DOI: https://dx.doi.org/10.3390/su151511525]

63. Esfandiari, M.; Jabari, S.; McGrath, H.; Coleman, D. Flood mapping using random forest and identifying the essential conditioning factors; a case study in Fredericton, New Brunswick, Canada. ISPRS Ann. Photogramm. Remote Sens. Spat. Inf. Sci.; 2020; 3, pp. 609-615. [DOI: https://dx.doi.org/10.5194/isprs-annals-V-3-2020-609-2020]

64. Chen, Z.; Sun, Z.; Zhang, H.; Zhang, H.; Qiu, H. Aboveground Forest Biomass Estimation Using Tent Mapping Atom Search Optimized Backpropagation Neural Network with Landsat 8 and Sentinel-1A Data. Remote Sens.; 2023; 15, 5653. [DOI: https://dx.doi.org/10.3390/rs15245653]

65. Speiser, J.L.; Miller, M.E.; Tooze, J.; Ip, E. A comparison of random forest variable selection methods for classification prediction modeling. Expert Syst. Appl.; 2019; 134, pp. 93-101. [DOI: https://dx.doi.org/10.1016/j.eswa.2019.05.028]

66. Mälicke, M. SciKit-GStat 1.0: A SciPy-flavored geostatistical variogram estimation toolbox written in Python. Geosci. Model Dev.; 2022; 15, pp. 2505-2532. [DOI: https://dx.doi.org/10.5194/gmd-15-2505-2022]

67. Musthafa, M.; Singh, G.; Kumar, P. Comparison of forest stand height interpolation of GEDI and ICESat-2 LiDAR measurements over tropical and sub-tropical forests in India. Environ. Monit. Assess.; 2023; 195, 71. [DOI: https://dx.doi.org/10.1007/s10661-022-10657-w]

68. Li, Y.; Li, M.; Liu, Z.; Li, C. Combining kriging interpolation to improve the accuracy of forest aboveground biomass estimation using remote sensing data. IEEE Access; 2020; 8, pp. 128124-128139. [DOI: https://dx.doi.org/10.1109/ACCESS.2020.3008686]

69. Byrd, K.B.; O’Connell, J.L.; Di Tommaso, S.; Kelly, M. Evaluation of sensor types and environmental controls on mapping biomass of coastal marsh emergent vegetation. Remote Sens. Environ.; 2014; 149, pp. 166-180. [DOI: https://dx.doi.org/10.1016/j.rse.2014.04.003]

70. Zhen, X.; Lu, L. Geostatistics (Modern Spatial Statistics); Science Press: Beijing, China, 2018.

71. Meiyan, L.; Sheng, N.; Cheng, W.; Xiaohuan, X.; Feng, C.; Baokun, F. Forest volume inversion based on ICESat-2 and Sentinel-2A data. Remote Sens. Nat. Resour.; 2024; 36, pp. 210-216. (In Chinese)

72. Meng, P.; Wang, H.; Qin, S.; Li, X.; Song, Z.; Wang, Y.; Yang, Y.; Gao, J. Health assessment of plantations based on LiDAR canopy spatial structure parameters. Int. J. Digit. Earth; 2022; 15, pp. 712-729. [DOI: https://dx.doi.org/10.1080/17538947.2022.2059114]

73. Meng, J.; Du, X.; Wu, B. Generation of high spatial and temporal resolution NDVI and its application in crop biomass estimation. Int. J. Digit. Earth; 2013; 6, pp. 203-218. [DOI: https://dx.doi.org/10.1080/17538947.2011.623189]

74. Liu, Y.; Chen, X.; Yang, J.; Li, L.; Wang, T. Snow avalanche susceptibility mapping from tree-based machine learning approaches in ungauged or poorly-gauged regions. Catena; 2023; 224, 106997. [DOI: https://dx.doi.org/10.1016/j.catena.2023.106997]

75. Bulut, S. Machine learning prediction of above-ground biomass in pure Calabrian pine (Pinus brutia Ten.) stands of the Mediterranean region, Türkiye. Ecol. Inf.; 2023; 74, 101951. [DOI: https://dx.doi.org/10.1016/j.ecoinf.2022.101951]

76. Zhang, L.; Zhang, X.; Shao, Z.; Jiang, W.; Gao, H. Integrating Sentinel-1 and 2 with LiDAR data to estimate aboveground biomass of subtropical forests in northeast Guangdong, China. Int. J. Digit. Earth; 2023; 16, pp. 158-182. [DOI: https://dx.doi.org/10.1080/17538947.2023.2165180]

77. Silveira, E.M.; Silva, S.H.G.; Acerbi-Junior, F.W.; Carvalho, M.C.; Carvalho, L.M.T.; Scolforo, J.R.S.; Wulder, M.A. Object-based random forest modelling of aboveground forest biomass outperforms a pixel-based approach in a heterogeneous and mountain tropical environment. Int. J. Appl. Earth Obs. Geoinf.; 2019; 78, pp. 175-188. [DOI: https://dx.doi.org/10.1016/j.jag.2019.02.004]

78. Chen, L.; Wang, Y.; Ren, C.; Zhang, B.; Wang, Z. Assessment of multi-wavelength SAR and multispectral instrument data for forest aboveground biomass mapping using random forest kriging. For. Ecol. Manag.; 2019; 447, pp. 12-25. [DOI: https://dx.doi.org/10.1016/j.foreco.2019.05.057]

79. Du, H.; Zhou, G.; Fan, W.; Ge, H.; Xu, X.; Shi, Y.; Fan, W. Spatial heterogeneity and carbon contribution of aboveground biomass of moso bamboo by using geostatistical theory. J. Plant Ecol.; 2010; 207, pp. 131-139. [DOI: https://dx.doi.org/10.1007/s11258-009-9659-3]

80. Lu, Y. Spatial-Temporal Co-Kriging Interpolation Method for Air Pollution Index Analysis; Chinese Academy of Surveying and Mapping: Beijing, China, 2018.

81. Zhou, Y.; Xie, B.; Li, M. Regional forest aboveground biomass mapping based on random forest and Kriging method—A case study of forest in northern Guangdong. J. Nanjing For. Univ. (Nat. Sci. Ed.); 2024; 48, 169.

82. Li, J.; Li, C.; Yin, Z. Kriging interpolation method based on ArcGIS and its application. Bull. Surv. Map.; 2013; 9, 87–90+97

83. Dowd, P.A.; Pardo-Igúzquiza, E. The many forms of co-kriging: A diversity of multivariate spatial estimators. Math. Geosci.; 2024; 56, pp. 387-413. [DOI: https://dx.doi.org/10.1007/s11004-023-10104-7]

84. Lu, Y.; Wang, L.; Qiu, A. A co-Kriging interpolation method based on principal component analysis. Bull. Surv. Map.; 2017; 11, 51.

85. High Carbon Stock Mapping at Large Scale with Optical Satellite Imagery and Spaceborne LIDAR. Available online: https://arxiv.org/abs/2107.07431 (accessed on 20 January 2024).

86. Dang, A.T.N.; Nandy, S.; Srinet, R.; Luong, N.V.; Ghosh, S.; Kumar, A.S. Forest aboveground biomass estimation using machine learning regression algorithm in Yok Don National Park, Vietnam. Ecol. Inf.; 2019; 50, pp. 24-32. [DOI: https://dx.doi.org/10.1016/j.ecoinf.2018.12.010]