Aboveground Forest Biomass Estimation Using Tent Mapping Atom Search Optimized Backpropagation Neural Network with Landsat 8 and Sentinel-1A Data

Abstract

Accurate forest biomass estimation serves as the foundation of forest management and holds critical significance for a comprehensive understanding of forest carbon storage and balance. This study aimed to integrate Landsat 8 OLI and Sentinel-1A SAR satellite image data and selected a portion of the Shanxia Experimental Forest in Jiangxi Province as the study area to establish a biomass estimation model by screening influencing factors. Firstly, we extracted spectral information, vegetation indices, principal component features, and texture features within 3 × 3-pixel neighborhoods from Landsat 8 OLI. Moreover, we incorporated Sentinel-1’s VV (vertical transmit–vertical receive) and VH (vertical transmit–horizontal receive) polarizations. We proposed an ensemble AGB (aboveground biomass) model based on a neural network. In addition to the neural network model, namely the tent mapping atom search optimized BP neural network (Tent_ASO_BP) model, partial least squares regression (PLSR), support vector machine (SVR), and random forest (RF) regression prediction techniques were also employed to establish the relationship between multisource remote sensing data and forest biomass. Optical variables (Landsat 8 OLI), SAR variables (Sentinel-1A), and their combinations were input into the four prediction models. The results indicate that Tent_ ASO_ BP model can better estimate forest biomass. Compared to pure optical or single microwave data, the Tent_ASO_BP model with the optimal combination of optical and microwave input features achieved the highest accuracy. Its R² was 0.74, root mean square error (RMSE) was 11.54 Mg/ha, and mean absolute error (MAE) was 9.06 Mg/ha. Following this, the RF model (R² = 0.54, RMSE = 21.33 Mg/ha, MAE = 17.35 Mg/ha), SVR (R² = 0.52, RMSE = 17.66 Mg/ha, MAE = 15.11 Mg/ha), and PLSR (R² = 0.50, RMSE = 16.52 Mg/ha, MAE = 12.15 Mg/ha) models were employed. In conclusion, the BP neural network model improved by tent mapping atom search optimization algorithm significantly enhanced the accuracy of AGB estimation in biomass studies. This will provide a new avenue for large-scale forest resource surveys.

1. Introduction

Aboveground biomass (AGB) is one of the important indicators for assessing the carbon sequestration capacity of forests [1]. The rapid and accurate estimation of forest biomass is a highly challenging research topic, holding significant implications for estimating forest carbon storage and sustainable development. The most accurate method for estimating forest biomass is based on field measurements, but collecting data through field measurements is both time-consuming and labor-intensive, especially in rugged terrains such as deep valleys or steep slopes [2]. Remote sensing systems, including optical and active sensors, have been proven to be effective alternative means for measuring and monitoring forest aboveground biomass at various scales [3,4]. By combining multisource remote sensing data with various models and sampling surveys, rapid data acquisition for forest resources, spatial structural distribution, and dynamic changes can be achieved, enabling quantitative inversion of forest measurement parameters [5,6]. Recent research has predominantly focused on combining optical sensors and synthetic aperture radar (SAR) for estimating and mapping forest AGB [7,8,9].

Landsat 8 (L8) and Sentinel-1 (S1) have been universally used for natural resource research, long a major source of data for AGB research. Launched in February 2013, Landsat 8 was designed by the National Aeronautics and Space Administration (NASA). With a multispectral operational land imager (OLI) sensor and thermal infrared sensor (TIRS), L8 provides data through 11 spectral bands with an image acquisition cycle of 16 days. The combination of different spectral information received by the sensor into vegetation indices provides strong data support for studying the spatial distribution of vegetation biomass in ecosystems. Based on L8 images, Ding et al. used the normalized difference vegetation index (NDVI), the wide dynamic range vegetation index (WDRVI), the chlorophyll index (CI), and Band6 and Band7 to build the regression kriging (RK) method to estimate grassland carbon stocks in Northeast China [10]. Shu et al. [11] used random forest regression (RFR), the nearest neighbor (K-NN) method, and partial least squares regression (PLSR) to study the aboveground biomass of white pine in Shangri-La City; the results indicate that the method based on the RFR biomass model worked best, and the relative accuracy reached up to 98.39%.

Although Landsat 8 (L8) imagery has shown promising results in the estimation of forest aboveground biomass (AGB), its capability to estimate higher biomass levels is limited as it is more sensitive to canopy density/coverage [12]. The issue of biomass saturation with low-to-medium spatial resolution optical sensors is a well-recognized problem [13,14,15], largely attributed to cloud cover and significant atmospheric aerosol interference [16]. On the contrary, synthetic aperture radar (SAR) sensors such as Sentinel-1, ALOS PALSAR, and TerraSAR-X have the capability to penetrate cloud cover and forest canopies, making them suitable for estimating forest AGB. Sentinel-1 consists of polar-orbiting satellites A and B, both equipped with SAR sensors, making them active microwave remote sensing satellites. The advantage of microwave remote sensing lies in its ability to overcome the influence of weather conditions and to collect data day and night, enabling continuous monitoring under all weather conditions; it contains vertical–vertical (VV) and vertical–horizontal (VH) polarizations. VH, in particular, has high sensitivity to the dynamics of vegetation density and structure [17]. Eduarda et al. [18] used indicators from Sentinel-1, along with latitude and longitude information, to build a model. The results indicate that VH polarization and longitude indicators are overall the most important predictors for AGB estimation.

The existing research indicates that the collaborative use of optical and SAR (synthetic aperture radar) data can provide more accurate estimates of various forest physical parameters compared to using a single sensor. More specifically, Zhang et al. [19] used Sentinel-1 and Sentinel-2 data to estimate forest AGB in northeastern Conghua, Guang-dong Province, China, through the implementation of the random forest regression (RFR) technique. The use of fused data yielded the highest accuracy, with a coefficient of determination (R²) = 0.72, compared to the individual analysis of optical and SAR data.

In addition, it is critical to select an appropriate algorithm to establish AGB estimation models. Linear regression methods have been the most commonly used approach in forest biomass studies in the past [20,21]. Classical statistical regression methods are simple and easy to implement. However, these classical statistical regression methods may not effectively capture the complex nonlinear relationships between forest aboveground biomass (AGB) and remote sensing data. To overcome this challenge, machine learning methods such as k-nearest neighbors (KNN), random forests (RFs), support vector machines (SVMs), and backpropagation neural networks (BPNNs) have been applied to research on forest AGB [22,23,24,25]. The KNN algorithm assigns a weighted average of the response values of the most similar neighbors to the target unit of interest, taking into account the similarity of measurement in feature space [26]. Henrik et al. [27] successfully mapped large areas of forest aboveground biomass and stem volume by combining TanDEM-X, Sentinel-2, and field reference data using the KNN algorithm. RF constructs multiple decision trees through bootstrapping and combines their prediction results to improve the accuracy and robustness of the model; the bootstrap is a statistical approach used to quantify the uncertainty associated with a given estimator [28]. The SVR method is aimed at estimating the correlation between a range of predictive variables and the target variable. The SVR has demonstrated its superiority in predicting vegetation parameters by relating spectral metrics to field samples [29]. Previous studies have shown that the artificial neural network algorithm based on backpropagation demonstrates excellent performance in biomass estimation modeling [30]. Zhou et al. [31] utilized high-resolution RGB imagery to calculate 17 vegetation indices and constructed linear, exponential, and backpropagation neural network models based on these indices to estimate AGB. The results indicated that the constructed BPNN model yielded the best estimation performance, with a validation R² of 0.68. Ying et al. utilized the atom search optimized BP neural network (ASO-BP) and the traditional backpropagation neural network (BP) to construct an aboveground biomass (AGB) inversion model for the Liupan Shan Forest Area. Through comparison, it was observed that the ASO-BP inversion model showed an average improvement of 6% in R², with a decrease of 2.4 Mg/ha in RMSE [32].

Combining multisource remote sensing data with forest inventory data for regional forest biomass research could reduce the uncertainty in the estimation, provides a more consistent spatial and temporal analyses [33,34]. In this study, through the integration of existing machine learning techniques, we combined forest field survey data with multisource remote sensing data. Additionally, an improved neural network algorithm was employed for the first time to estimate forest biomass, aiming to enhance the accuracy of biomass estimation at the regional scale. Our objectives were as follows: (1) Test the reliability from field-based AGB and the multisource remote sensing parameter as training and validation samples for subsequent satellite-based modeling. (2) Explore the impact of an improved tent mapping atom search optimized BP neural network algorithm, compared to traditional methods, including partial least squares regression (PLSR), random forest regression (RFR), and support vector machine regression (SVR), on the accuracy of AGB prediction, and evaluate the performance metrics using the root mean square error (RMSE), mean absolute error (MAE), and coefficient of determination (R²). (3) Combine parameter and nonparameter methods to select variables, establish the optimal biomass model for the study area’s forest, and generate biomass maps.

2. Materials and Methods

2.1. Study Area

The study area, the Subtropical Forestry Experimental Center, Chinese Academy of Forestry, is located in the central Jiangxi Province in Southwestern China (Figure 1). The geographical coordinates are 114°37′55″–114°43′46″ E, 27°43′9″–27°45′25″ N. The total area is 1861.3 hm², and it has a forest coverage rate of 95.6% [35,36]. This area has a typical south subtropical humid continental monsoon climate. The altitude is 220–1092 m, the terrain has large ups and downs, the annual average temperature is 17.2 °C, and the annual average precipitation is 1950 mm. Rainfall is mostly concentrated in March to June. The average canopy density of trees in the study area is above 0.7, and there are many tree species, including Cunninghamia lanceolata, Schima superba, Pinus massoniana, Liriodendron chinense, Liquidambar formosana, and Machilus pauhoi Kaneh.

2.2. Field Data

The field data were collected at the end of November 2019 (Figure 2) from representative small plots within the forest area, with an average stand age of 30 years. During this season, the growth traits of the trees were in a stable phase [37,38]. A total of 4094 individual tree measurements were collected, covering an area of 75.52 acres. The collected data for each plot include tree species information, individual latitude and longitude coordinates, elevation, tree species category, tree height, diameter at breast height (DBH), and crown width. A 5 cm diameter was used as the starting measurement size, and a DBH tape measure was employed to measure the DBH of each individual tree, with the precision maintained within 1 to 2 cm. Tree height and the east–west and north–south crown widths were measured using a laser rangefinder and a measuring tape. To ensure the accuracy of individual tree coordinates, Real-time Kinematic (RTK) carrier phase differential technology was used for positioning and recording the coordinates of each individual tree.

The calculation of individual tree biomass within the plots was carried out using China’s primary forest biomass modeling manual [39]. Given the spatial resolution of Landsat 8 OLI optical remote sensing at 30 m × 30 m, individual tree positions were geolocated onto the remote sensing imagery. This facilitated the computation of single tree biomass values for different tree species within each pixel area. The biomass of each individual tree species was then summed up to obtain the corresponding aboveground biomass (AGB) at the plot scale, denoted as $W_p$; considering that tree growth conditions can vary depending on different regions and climatic conditions, biomass calculation equations for various dominant tree species were utilized in conjunction with the location of the forest stand. These equations are presented in

Table 1. Aboveground biomass model of main tree species, where W represents the biomass to be calculated, measured in kilograms (kg); DBH stands for the individual tree’s diameter at breast height, measured in centimeters (cm); and H represents the height of the individual tree, measured in meters (m).

Species of Trees	Aboveground Biomass Model	Reference
Cunninghamia lanceolata	$W_1=0.0183{\left(DB{H}^2\times H\right)}^{0.9663}$	[40]
Schima superba	$W_2=0.1204\times DB{H}^{2.0644}\times H^{0.3826}$	[41]
Pinus massoniana	$W_3=0.0973{\left(DB{H}^2\times H\right)}^{0.8284}$	[42]
Liriodendron chinense	$W_4=0.0253{\left(DB{H}^2\times H\right)}^{0.9747}$	[43]
Liquidambar formosana	$W_5=0.0252{\left(DB{H}^2\times H\right)}^{0.9614}$	[44]
Machilus pauhoi Kaneh	$In\left(W_6\right)=-3.0557+0.9429\times In\left(DB{H}^2\times H\right)$	[45]
Plotted biomass calculation	$W_p=\sum \left(W_1+W_2+\dots +W_n\right)$	[39]

.

Combining remote sensing-scale clipping processing, a total of 50 sample plots of biomass were established, with plot biomass ranging from 14.92 Mg/ha to 126.96 Mg/ha, covering forest areas from sparse plots to dense plots. Among them, tree density represents the tree density converted to trees per hectare based on the number of trees corresponding to each 30 m × 30 m plot. The specific survey information is shown in

Table 2. Field sample survey dataset (different plots and biomass).

Forest Indicators	Field Sample (n = 50)
	Max	Min	Mean
Biomass (Mg/ha)	126.96	14.92	78.08
Mean H (m)	17.29	8.18	12.53
Mean DBH (cm)	25.11	14.02	18.6
Tree density (no. trees/ha)	1244.44	111.11	711.11

.

After estimating the forest biomass, the biomass–carbon conversion factor was used to transform forest biomass into the corresponding forest carbon storage. The formula is as follows:

C = W × CF

(1)

where C denotes the forest carbon storage, W denotes the forest biomass, and CF denotes the conversion factor. Research has shown that this conversion factor is associated with tree species, but it remains relatively stable within a certain range (0.47–0.55) [46]. In this study, the median conversion factor (0.5) was chosen to offset the positive and negative errors in carbon stock estimation [47].

2.3. Landsat 8 Data and Preprocessing

Landsat 8 OLI data have been widely used in biomass estimation studies [10]. The Landsat 8 multispectral images in this study were downloaded from USGS Earth Explorer (https://earthexplorer.usgs.gov/, accessed on 4 May 2023), covering bands B2–B7. In consideration of synchronizing the acquisition time of remote sensing imagery with survey data, a scene from 14 November 2019, with clear sky conditions and a path/row number of 122/41 was selected.

In ENVI 5.3 software, Landsat 8 OLI remote sensing images undergo radiometric calibration and atmospheric correction using the FLAASH model. Radiometric calibration is a process that converts Digital Numbers (DNs) into absolute radiance values, thus eliminating errors inherent to the sensor itself. Atmospheric correction is applied to mitigate the effects of aerosol scattering, allowing for the retrieval of genuine model parameters such as surface reflectance and land surface temperature. Additionally, vector clipping is performed in conjunction with the spatial coordinates of the Shanxia Experimental Forest to refine the data.

2.4. Synthetic Aperture Radar (SAR) Data and Preprocessing

Sentinel-1 is a C-band SAR remote sensing satellite that launched into orbit in 2014. The main objective is to provide all-weather, all-season, high-resolution Earth surface observation data. There are four imaging models (Stripmap (SM), Interferometric Wide swath (IW), Extra Wide swath (EW), and Wave (WV)), along with different polarization modes, including single-polarization (HH, VV, HV, and VH) and multipolarization (HH&HV, VV&VH).

The SAR data used in this study were Ground Range Detected (GRD) products with a spatial resolution of 5 m × 20 m at multiple viewing angles, with the VV (vertical transmit–vertical receive) and VH (vertical transmit–horizontal receive) polarizations, obtained on 16 November 2019. The dataset was preprocessed using the Sentinel Applications Platform (SNAP 9.0.0) of the ESA software, orbit correction, thermal noise removal, radiometric calibration, speckle filtering, and terrain correction (resampled using SRTM 1Sec HGT) for Sentinel-1 SAR data. In order to mitigate the influence of low values in the images, we converted the preprocessing results to dB units (Equation (2)). Subsequently, all the S-1 data were resampled to 30 m resolution using the bilinear interpolation method and resized based on the spatial coverage of the study site.

dB = 10 × log₁₀(N)

(2)

where N is the value per pixel from the preprocessed SAR images.

3. Methods and Modelling

Figure 3 presents an overview of the AGB modeling process; it can be divided into the following three steps.

The first step is data preprocessing and parameter extraction. In the second step, different screening rules are applied to filter parametric and nonparametric predictors, and 4 different sets of experimental plans are designed based on different data source characteristics. The third step is to build a model; for the four sets of experimental plans and field sampling datasets in the second stage, we randomly selected 70% of the data as the training dataset, built models based on PLSR, SVR, RF, and Tent_ASO_BP algorithms, and compared the R², RMSE (Mg/ha), and MAE (Mg/ha), thereby determining the optimal model. The optimal model, along with the selected explanatory variables, was then used to generate an aboveground biomass (AGB) map of the study area.

3.1. The Extraction of Feature Variables from Remote Sensing Data

Three classes of optical variables were calculated from the Landsat 8 OLI reflectance images: spectral bands, spectral indices, and textural features. The extracted spectral bands included blue, green, red, near-infrared, SWIRI, and SWIRII. Based on the relevant literature on the sensitivity of optical features to biomass, 11 spectral vegetation indices were calculated, including the normalized difference vegetation index (NDVI), green ratio vegetation index (GRVI), difference vegetation index (DVI), enhanced vegetation index (EVI), simple ratio (SR), soil-adjusted vegetation index (SAVI), optimized soil-adjusted vegetation index (OSAVI), specific leaf area vegetation index (SLAVI), albedo, atmospherically resistant vegetation index (ARVI), and normalized difference infrared index (NDII). The principal component analysis involved dimensionality reduction and feature extraction of the base bands. The tasseled cap transformation (TCT) is a widely used linear transformation commonly employed for characterizing forest structural conditions [48,49]. In this biomass study, three transformations, namely tasseled cap brightness (TCB), tasseled cap greenness (TCG), and tasseled cap wetness (TCW), were utilized.

In addition, we extracted texture features from the image using the Gray-Level Co-occurrence Matrix (GLCM). Previous research has indicated that the best AGB estimation accuracy can be achieved using a 3 × 3 window size [50]. Principal component analysis (PCA) conducted on the Landsat 8 OLI image bands revealed that Band 2 contains rich feature information. As a result, eight texture feature values were extracted from Band 2, including mean (ME), variance (VAR), contrast (CON), entropy (ENT), dissimilarity (DIS), second moment (SM), homogeneity (HO), and correlation (COR).

Taking into account the dual-polarization characteristics of Sentinel-1A data, and following previously published research [51,52], we chose to use two polarization bands (VV, VH), their ratio (VH/VV), their difference (VH-VV), and their sum (VH+VV) as inputs for SAR.

Table 3. Landsat 8 OLI and Sentinel-1A spectral and textural features used in the modeling ($L$ is 0.5; $r$ is 1).

Source	Type	Variable	Full Name
Landsat 8	Spectral bands	Blue (B2)
		Green (B3)
		Red (B4)
		Near-infrared (B5)
		SWIRI (B6)
		SWIRII (B7)
	Spectral indices	NDVI	Normalized difference vegetation index, $\frac{B_5-B_4}{B_5+B_4}$
		GRVI	Greenness ratio vegetation index, $\frac{B_5}{B_3}$
		DVI	Difference vegetation index,$B_5-B_4$
		EVI	Enhanced vegetation index,$\frac{2.5\times \left(B_5-B_4\right)}{B_5+6B_4-7.5B_2+1}$
		SR	Simple ratio vegetation index, $\frac{B_5}{B_4}$
		SAVI	Soil-adjusted vegetation index, $\frac{(B_5-B_4)}{B_5+B_4+L}\left(1+L\right)$
		OSAVI	Optimized soil regulation vegetation index, $\frac{B_5-B_4}{B_5+B_4+L}$
		SLAVI	Specific leaf area vegetation index, $\frac{B_5}{B_4}+B_7$
		Albedo	Surface reflectance,$B_2+B_3+B_4+B_5+B_6+B_7$
		ARVI	Atmospherically resistant vegetation index, $\frac{B_5-\left(B_4-\gamma \left(B_2-B_4\right)\right)}{B_5+\left(B_4-\gamma \left(B_2-B_4\right)\right)}$
		NDII	Normalized difference infrared index,$\frac{B_5-B_6}{B_5+B_6}$
		PCA	Pca1, Pca2, Pca3, Pca4, Pca5, Pca6, Pca7
		TCB	Tasseled cap brightness
		TCG	Tasseled cap greenness
		TCW	Tasseled cap wetness
	Textural features	ME	Mean,${\sum }_{i,j=0}^{N-1}i{P}_{i,j}$
		VAR	Variance, ${\sum }_{i,j=0}^{N-1}{P}_{i,j}\left(1-{\mu }_i\right)$
		HO	Homogeneity, ${\sum }_{i,j=0}^{N-1}i\frac{P_{i,j}}{1+{\left(i-j\right)}^2}$
		CON	Contrast, ${\sum }_{i,j=0}^{N-1}i{P}_{i,j}{\left(i-j\right)}^2$
		DIS	Dissimilarity, ${\sum }_{i,j=0}^{N-1}i{P}_{i,j}\left|i-j\right|$
		ENT	Entropy, ${\sum }_{i,j=0}^{N-1}i{P}_{i,j}\left(-In{P}_{i,j}\right)$
		SM	Second moment, ${\sum }_{i,j=0}^{N-1}i{P}_{i,j}{}^2$
		COR	Correlation, ${\sum }_{i,j=0}^{N-1}i\frac{{\sum }_{i,j=0}^{N-1}i,j{P}_{i,j}-{\mu }_i{\mu }_j}{{\sigma }_i^2{\sigma }_j^2}$
		${\mu }_i={\displaystyle \sum _{i=0}^{N-1}}i{\displaystyle \sum _{j=0}^{N-1}}{P}_{i,j}; {\mu }_j={\displaystyle \sum _{j=0}^{N-1}}j{\displaystyle \sum _{j=0}^{N-1}}{P}_{i,j}; {\sigma }_i^2={\displaystyle \sum _{i=0}^{N-1}}{\left(i-{\mu }_i\right)}^2{\displaystyle \sum _{j=0}^{N-1}}{P}_{i,j};$ ${\sigma }_j^2={\displaystyle \sum _{j=0}^{N-1}}{\left(j-{\mu }_j\right)}^2{\displaystyle \sum _{j=0}^{N-1}}{P}_{i,j}$
Sentinel-1A		VV	Vertical transmit–vertical channel
		VH	Vertical transmit–horizontal channel
		VH/VV	Cross-polarized ratio
		VH-VV	Polarization difference
		VH+VV	Polarization sum

shows a list of all the variables (metrics extracted from L8 and S-1 data) employed in the study.

3.2. Important Features Identified for Forest AGB Prediction

Feature selection has a significant impact on the performance of regression models by eliminating irrelevant variables and addressing multicollinearity [53]. Proper feature combinations can notably enhance the estimation accuracy of an AGB model. The simple Spearman correlation, which assesses monotonic relationships, linear or nonlinear, has sometimes been used, but it has specific limitations in nonlinear relationship scenarios [54]. In this study, two feature selection methods are provided. The parametric model selects feature variables based on the variable projection importance in partial least squares regression (PLSR), with all $V_{IP}$ (Variable Importance in Projection) thresholds set at 1. The nonparametric model employs the feature importance measure of random forests to select feature variables. All variables are included in the model, and the relative importance of variables is calculated [55]. Variables with a significant increase in root mean square error percentage increase (%IncMSE) are chosen. The selected variable combinations are then used for further modeling.

Feature variable importance calculations were carried out using both SIMCA 14.1 software and the randomForest package in R version 4.2.3.

3.3. Regression Algorithms

In our study, four different regression methods were employed to map AGB, including both parametric and nonparametric approaches. These methods are PLSR, SVR, RF, and Tent_ASO_BP. Specifically, SVR and Tent_ASO_BP optimization models were implemented using MATLAB R2022a, while the RF model was developed using the R programming language (R 4.2.3). The description of each algorithm is provided below.

3.3.1. Partial Least Squares Regression

Partial least squares regression (PLSR), as a parameterized regression model, aims to predict a linear regression model between independent and dependent variables. This prediction is achieved by extracting a set of orthogonal factors with the best predictive capability from the independent variables to meet modeling requirements. PLSR, as an extension of multivariate linear regression models, builds upon the foundation of ordinary multivariate regression analysis by incorporating the concepts of principal component analysis and canonical correlation analysis. It uses cross-validation techniques to assess the quality and effectiveness of its regression model [56]. Furthermore, PLSR includes an important variable selection process known as Variable Importance in Projection ($V_{IP}$) (Equation (3)), which is used to determine the explanatory power of individual independent variables for the dependent variable.

$$V_{IP}=\sqrt{\frac{k}{{{\displaystyle \sum }}_{h=1}^n{r}^2\left(y,c_h\right){{\displaystyle \sum }}_{h=1}^n{r}^2\left(y,c_h\right)w_{hj}^2}}$$

(3)

where $k$ represents the number of independent variables, $c_h$ denotes the principal component extracted from the independent variables, $r\left(y,c_h\right)$ stands for the correlation coefficient between the dependent variable and the principal component, and $w_{hj}$ represents the weight of the independent variable in the principal component.

The larger the importance of an independent variable for the dependent variable, the larger its $V_{IP}$ value, and vice versa [57]. Researchers have suggested using a threshold value of 0.8 to distinguish between important and unimportant variables [58]. In this study, parameterized regression modeling selected projection importance as the variable selection criterion.

3.3.2. Support Vactor Regression

According to statistical learning theory, support vector regression (SVR) is one of the most effective machine learning techniques [59,60], and it has demonstrated superiority in its generalization ability [23]. As a type of nonparametric model, it can be employed in a manner similar to supervised learning to address regression problems [61]. The principle behind it is to use a nonlinear mapping kernel function to map input vectors into a high-dimensional feature space to achieve regression [62]. The commonly used kernel functions include the polynomial kernel function, radial basis function (RBF), and sigmoid kernel function. In this study, the radial basis kernel function was used, which has been frequently employed in previous research for optimizing forest aboveground biomass (AGB). It has fewer parameters, with the model primarily determined by two hyperparameters: the penalty coefficient C and the gamma parameter value g. The penalty coefficient C represents the model’s tolerance to errors, and both excessively large and small C values can lead to a decrease in the model’s generalization ability. The gamma parameter, implicitly, determines how the data are mapped to a new feature space, and its magnitude directly affects the speed of model training and prediction. Therefore, a grid search and five-fold cross-validation parameter optimization method was employed to search for the most suitable combinations of C and g parameters within a predefined threshold range, maximizing the R² to ensure optimal model accuracy. Through experimentation, it was found that for the Landsat 8 dataset, the optimal SVR model parameters were C = 0.5 and g = 0.03, while for the Sentinel-1A dataset, the optimal parameters were C = 8 and g = 0.015. For the S-L8-S1 dataset, the optimal values were found to be C = 8 and g = 0.031.

3.3.3. Random Forest Regression

The RF algorithm is considered to be the most commonly used and effective algorithm. It belongs to the ensemble learning methods, roughly divided into bagging and boosting methods [62]. For regression, RF employs the boosting method to randomly select a subset of features with replacement from the sample set and constructs decision trees based on this subset of features. These decision trees are combined to form a forest, and the results of all the decision trees are weighted and averaged to produce the final prediction for the sample. During the construction of decision trees, it is necessary to rank the importance of the features. The importance of features in random forests is measured by two indicators. The first indicator assesses the impact of independent variables on the purity of tree nodes using out-of-bag data. This is calculated by the sum of squared residuals, with higher values indicating greater importance. The second metric involves differences based on the mean squared error; comparing the fitting results of out-of-bag data before and after shuffling into decision trees, the more important features are determined. In the feature selection process, such feature data are input into the random forest model as independent variables. Within the random forest regression algorithm, two important parameters exist: the number of trees (ntree) and the number of features at each split point (mtry). Among them, ntree controls the number of trees that adjust the regression error, while mtry determines the number of features used for node splitting. Only the coordination of the above hyperparameters can train a well-performed model [63]. Random forests are commonly used for modeling with single predictive variables and combinations of multiple predictive variables.

3.3.4. Tent Mapping Atom Search Optimization BP Neural Network

The BP neural network is a type of multilayer feedforward neural network based on error backpropagation. It consists of an input layer, hidden layer(s), and an output layer, with both the input and output layers having only one layer, while there can be one or more hidden layers. The basic BP neural network algorithm comprises two stages: signal forward propagation and error backward propagation. In the first stage, input data are passed through the network, calculated layer by layer, and transmitted to the output layer to generate the network’s prediction. This process involves input data, including weight and bias parameters. The second stage involves the error being propagated backward from the output layer, passing through the hidden layer(s), and reaching the input layer. During this process, adjustments are made sequentially to the weights and biases of each layer to reduce prediction errors. This is known as the backpropagation of errors. Throughout the entire training process, forward propagation and backward propagation alternate until the neural network’s performance meets the desired criteria or reaches the training stop condition.

The atom search optimized algorithm is a novel physics-inspired metaheuristic optimization algorithm that was developed in 2018. In this algorithm, the physical interpretation is that the position of each atom within the search space represents a solution measured by its mass. Based on the distances between atoms and the interaction forces generated by the Lennard–Jones potential, all atoms in the population either attract or repel each other. Lighter atoms are attracted towards heavier ones, and the optimization process involves updating the positions of these atoms throughout the entire search space to seek the optimal solution [64]. The basic steps of the algorithm include initializing the atom population, calculating fitness, computing gradients, updating atom positions, managing atom interactions, constraining positions, and iteratively optimizing until the algorithm outputs the best-found solution along with its corresponding fitness value. It is important to note that the acceleration of atoms arises from two components. One is the interaction force caused by the L-J potential (Equation (4)), which is in fact the vector sum of the attraction and the repulsion exerted from other atoms. Another one is the constraint force caused by the bond length potential, where the bond length potential is the weighted positional difference between each atom and the optimal atom [64].

Tent chaotic mapping is a nonlinear mapping method characterized by its ability to generate highly random sequences of numbers. Therefore, it is commonly used for initializing the weights and bias parameters of neural networks, effectively preventing the network from falling into local optima.

$$E\left(r\right)=4ϵ\left[{\left(\frac{\sigma }{\gamma }\right)}^{12}-{\left(\frac{\sigma }{\gamma }\right)}^6\right]$$

(4)

In the equation, $E\left(r\right)$ represents the interaction energy between two atoms, $\gamma$ is the Euclidean distance between the two atoms, $ϵ$ denotes the strength of attraction between molecules, and $\sigma$ represents the minimum distance between molecules.

BP neural networks can effectively handle nonlinear relationships among data [65]. However, due to their simulation process being somewhat similar to a black-box operation, their underlying mechanisms are difficult to describe for improving the accuracy and fit of the neural network model. In order to enhance the accuracy of model inversion and fitting, the tent mapping atom search optimized algorithm is introduced to optimize the BP neural network for forest aboveground biomass inversion [66]. This algorithm combines chaos theory and the atom search concept, generating a set of initial atom positions through tent chaotic mapping. These positions correspond to different points in the possible solution space. Next, atom search optimization is introduced into the BP network model, initializing the weights and thresholds of BP and calculating the fitness of the initial atoms. At this stage, each chaotic atom contains information such as weights, bias parameters, and learning rates. Subsequently, the atom search optimization iteration begins, continuously optimizing based on the speed generated by atom interactions and model error fitness, ultimately achieving a more precise and rapid discovery of globally optimal weight and threshold parameter solutions. The Tent_ASO_BP network structure is illustrated in Figure 4.

3.4. Experimental Design and Predictors

We used a total of 40 variables to estimate forest aboveground biomass (AGB) and designed four sets of experiments for the study (

Table 4. Experimental design for AGB estimation.

Experiment	Abbreviation	Description
1	L8	All spectral bands, vegetation indices, band transformations, and texture features.
2	S1	Polarization bands and corresponding sum, difference, and quotient bands.
3	L8-S1	All obtained L8 and S1 predictors.
4	S-L8-S1	Represents the predictions selected from the L8-S1 dataset using the different filtering rules in Figure 5.

). Among the selected predictors, (1) L8 included the following: 6 spectral bands, 11 vegetation indices, 7 principal component analysis features, 3 wavelet transforms, and 8 texture variables; (2) S1 included the following: 2 SAR polarization channels and 3 polarization indices. The extracted biomass values were appended to the calculated backscatter and optical information [19].

3.5. Model Assessment

To verify the estimation effect of the model, 70% (n = 35) of the fieldwork samples were randomly selected as training samples to build the model, and the remaining 30% (n = 15) were used for validation.

We used the root mean square error (RMSE, Mg/ha), coefficient of determination (R²), and mean absolute error (MAE, Mg/ha) to compare the performance of the selected machine learning techniques in forest AGB estimation. These metrics are widely employed in forest AGB modeling to assess the disparities between observed data and predicted forest AGB [67,68]. RMSE and MAE reflect the accuracy of the models, while R² is used to evaluate the goodness of fit; a higher R² indicates a better model [69].

The calculation methods for these statistical data are as follows:

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

(5)

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

(6)

$$MAE=\frac{1}{n}{\displaystyle \sum }_{i=1}^n\left|y_i-{\widehat{y}}_i\right|$$

(7)

where $y_i$ represents the measured biomass of the ith plot, ${\widehat{y}}_i$ represents the predicted biomass of the ith plot, and ${\overline{y}}_i$ represents the average biomass of n plots.

4. Results

4.1. Variable Selection Result

Using variable importance projection for feature selection, the results (Figure 5a) show significant correlations between various vegetation indices in optical remote sensing and biomass. SR exhibited the most pronounced correlation, with the highest explanatory power for biomass parameters, followed by NDII. In the radar bands, coefficients for both polarizations, VV (vertical–vertical) and VH (vertical–horizontal), as well as the sum of coefficients, demonstrated good explanatory power for AGB.

In the nonparametric feature selection, variable importance measures %IncMSE and the probability values of mean squared error increase (p-value) were computed. By default, the variables with p-values less than the significance level (typically 0.05) were considered as important predictive factors. As shown in Figure 5b, NDII emerged as the most critical variable, and we found that the third principal component factor (PCA3) contained more biomass-related information. Backscatter coefficient VH was found to be more suitable for AGB inversion compared to VV [70]. Features such as ARVI, b7, CON, DIS, DVI, HO, NDII, NDVI, OSAVI, PCA3, SAVI, SR, VH, and VH/VV contributed significantly, and will thus be used as input features for subsequent Tent_ASO_BP model development and AGB prediction.

The optimally selected features from parameterized and nonparameterized filtering were used as the original feature data for the PLSR model and the nonparametric model in the S-L8-S1 scheme of Section 4.2.

4.2. Comparison of AGB Experimental Models

Table 5. Modeling performance for diverse prediction algorithms and data sources.

Data	Algorithm	Features	R2	RMSE(Mg/ha)	MAE(Mg/ha)
L8	PLSR	35	0.42	23.99	21.65
	SVR	35	0.36	23.97	19.99
	RF	35	0.47	25.04	20.30
	Tent_ASO_BP	35	0.50	24.61	21.44
S1	PLSR	5	0.37	26.29	20.99
	SVR	5	0.29	25.56	20.47
	RF	5	0.33	21.54	17.42
	Tent_ASO_BP	5	0.48	25.12	20.73
L8-S1	PLSR	40	0.46	25.53	23.43
	SVR	40	0.45	16.44	13.93
	RF	40	0.51	20.32	16.31
	Tent_ASO_BP	40	0.60	14.34	11.23
S-L8-S1	PLSR	23	0.50	16.52	12.15
	SVR	14	0.52	17.66	15.11
	RF	14	0.54	21.33	17.35
	Tent_ASO_BP	14	0.74	11.54	9.06

Note: the best performing model (i.e., Tent_ASO_BP, with S-L8-S1 dataset) was put in bold.

displays the test results for 16 experiments conducted using L8 and S1 data. Before each experiment, the data sequence was randomly shuffled using the matlab randperm function to ensure the independence of the subsequent 70% training and 30% testing set splits (N = 15) and the model’s generalization performance. Among the four modeling schemes designed based on different feature variables, the Tent_ASO_BP model consistently outperformed the other models in terms of relative precision. The RF model followed, while the SVR model exhibited the poorest predictive performance. The performance of the parameterized PLSR model varied across different schemes. For instance, in the first three schemes, PLSR performed better overall than the SVR model. However, it showed relatively poorer performance in the modeling based on selected features, which is partly attributed to the feature selection rules. The RF model demonstrated relatively good performance thanks to the advantages of ensemble learning, which reduced the risk of overfitting by constructing multiple decision trees. Experiment 3, compared to experiments 1 and 2, indicated that combining optical remote sensing and radar data sources can effectively address the issue of data saturation associated with a single data source. It provided better estimation accuracy (an R² improvement of 0.10–0.18) and lower prediction error (the RMSE reduced by 1.22–10.78%). Among all models, Tent_ASO_BP had the best fitting performance, with S-L8-S1, based on combined data, achieving higher accuracy. According to the scatterplot of predicted AGB values (Figure 6), there are noticeable differences in modeling performance among the four prediction algorithms. Compared to other models, the Tent_ASO_BP model exhibited higher consistency between estimated values and ground observations. Furthermore, regardless of the prediction algorithm used, better estimation results were obtained compared to using optical or microwave predictions alone. The results indicate that the Tent_ASO_BP prediction algorithm, using feature-selected variables that combine optical and SAR data, achieved the highest R² and lowest RMSE and MAE values (0.74, 11.54 Mg/ha, 9.06 Mg/ha, respectively).

4.3. Forest AGB Predictive Mapping and Carbon Storage

Figure 7 displays the spatial distribution of AGB across the entire study area, utilizing optical and microwave features selected by nonparametric methods as inputs for the Tent_ASO_BP model. The AGB distribution ranged from 8.14 Mg/ha to 134.64 Mg/ha, with an average of 42.68 Mg/ha. The northern part of the study area was primarily characterized by the presence of lakes and rivers, which, due to the abundant presence of algae and organic matter, exhibit a deep green color. This affects the absorption and scattering properties of light, thereby influencing the signals received by the sensors, resulting in elevated values of NDVI, ARVI, and SR. According to the AGB map, the highest AGB values in the northern region of the study area fell in the range of 74 Mg/ha to 134 Mg/ha, indicating a phenomenon of biomass overestimation. In contrast, the majority of the study area (75%) showed lower AGB values of less than 50 Mg/ha, with an average AGB estimated at 37.37 Mg/ha. Furthermore, the southern side of the ridge, receiving abundant sunlight, exhibited a higher level of biomass compared to the northern side. Excluding a portion of the northern lake area, the biomass distribution ranged from 8.63 Mg/ha to 103.23 Mg/ha, with an average of 44.58 Mg/ha. Combining the AGB map generated by the Tent_ASO_BP model, the estimated biomass in the study area is 8.09 × 10⁴ t and there is a carbon storage level of 4.05 × 10⁴ t in the forested regions.

5. Discussion

5.1. The Contribution of L8 and S1 Derivatives in the Estimation of AGB

The results of variable feature selection (Figure 5) indicate that derivatives of L8 and S1 data are more suitable for predicting AGB than the original information, including backscatter and spectral bands. Similar findings have been reported in previous studies that estimated AGB using SAR and optical imagery [15,71,72]. For instance, Li et al. [71] found that backscatter indices (S1), spectral indices, and texture features (L8) outperformed multispectral (L8) and backscatter (S1) bands in AGB mapping. Liu et al. [72], in their investigation of biomass burning and biomass change before and after forest fires using multiple remote sensing sources, discovered that a combined model incorporating Landsat 8 data, Sentinel-1 data, and terrain variables achieved higher accuracy compared to models using the original bands (L8 and S1) and vegetation indices. Notably, they identified L8’s normalized burn ratio index (NBR) and normalized difference moisture index (NDMI) as the most critical predictors. Referencing Figure 5, in the case of nonparametric models, the selected indices of the quotient (VH/VV) provide more valuable information than the raw backscatter bands (VV). Among the various combinations of features from L8, the vegetation index NDII proved to be the optimal predictor for AGB estimation, followed by PCA3 and texture features (dissimilarity and homogeneity). This finding bears some resemblance to the conclusions drawn by Davide et al. [73], who studied mowing monitoring algorithms using six vegetation indices, including EVI, NDVI, and NDII. They found that NDII exhibited the best performance in monitoring when changes in ground biomass occurred due to mowing, achieving an overall accuracy of up to 93%. This suggests that the combination of SWIR and NIR bands plays a significant role in forest biomass prediction.

The results of this study indicate that, while Sentinel-1A data can mitigate the issue of optical image saturation caused by dense forest canopies, its relatively lower penetration capability in the C-band alone cannot serve as an effective predictor [71]. However, the combination of data yields better results compared to single remote sensing datasets. Simultaneously, it also demonstrates the feasibility of integrating Landsat 8 and Sentinel-1A, two crucial global freely available remote sensing observational datasets, in forest aboveground biomass (AGB) estimation.

5.2. The Performances of Different Prediction Models with Different Data Inputs

Neural networks belong to a class of powerful modeling algorithms capable of computing and predicting data, and they are more versatile than regression algorithms [74]. Relevant research employs backpropagation neural networks (BP) to estimate forest biomass in different regions [75,76,77] and confirms their superior performance in handling complex nonlinear relationships, with robustness and fault tolerance. In our study, compared to the other three methods, the Tent_ASO_BP model proved to be more efficient and superior. This can be attributed to Tent_ASO_BP’s strong adaptive adjustment strategy during global search. By mapping input variables into the network model and using tent mapping for initializing weights and bias parameters globally, the optimal parameter combination is obtained through iterative training, accelerating the convergence speed of the neural network and enabling accurate AGB predictions. However, there are certain limitations, as the atom search and tent mapping in the model can lead to increased computational complexity, especially in high-dimensional parameter spaces, which may require higher computational resources. In addition to Tent_ASO_BP, the overall predictive performance of the RF model is better than other prediction models. Other studies have also found that RF predictions exhibit good predictive accuracy [15,77,78]. This can be attributed to RF’s powerful ensemble learning capability and the stability of decision trees, which effectively reduce the risk of overfitting through the voting mechanism of multiple decision trees, achieving robustness against outliers. However, the RF method is less suitable for modeling with small sample sizes [79]. A possible explanation for the relatively inferior performance of support vector regression (SVR) compared to the RF and Tent_ASO_BP models in this study is the small sample size of the modeling data, and the influence of noisy variables on the construction of the optimal hyperplane, resulting in reduced model accuracy. Partial least squares regression (PLSR) combines the advantages of principal component analysis (PCA) and multivariate linear regression, and it can better handle high-dimensional data and multicollinearity issues. In single-data-source studies, the PLSR model demonstrates better accuracy than traditional machine learning methods. Although there is some improvement in performance after feature selection with parameterized rules, and a reduction in data dimensionality, the overall accuracy improvement in the PLSR model remains lower than that of machine learning methods. In future modeling, conducting a univariate analysis of AGB correlations, employing leave-one-out cross-validation for results on both the training and test sets, and experimenting with different hyperparameters to enhance the performance of AGB modeling may be considered.

In our study, we made the first attempt to introduce an improved backpropagation (BP) neural network model based on the tent mapping and atom search algorithm for estimating forest biomass. We compared this model with both parametric and nonparametric models. The results indicate a significant enhancement in the accuracy of the improved neural network model. Meanwhile, it is our belief that the BP neural network improved based on tent mapping and ASO atom optimization has universality. Previous studies have applied the Tent_ASO_BP algorithm in indoor positioning for smartphones [80], lithium battery health prediction [81], and construction safety management in the construction industry [82], all achieving favorable results. In the future, more optimization algorithms can be further improved and applied to research in the direction of forest biomass.

6. Conclusions

In this study, we explored the potential of combining SAR data and optical remote sensing data for estimating AGB. Based on our modeling approach, we draw the following conclusions: (1) The application of the improved backpropagation neural network learning model (Tent_ASO_BP) in estimating forest aboveground biomass (AGB) was explored for the first time. It produced the best estimation results compared to the other tested prediction techniques, including PLSR, SVR, and RF. (2) By comparing with the original dataset, it was found that variables selected based on two different screening rules provided strong data support for the optimal model. (3) Compared to optical or microwave data, the synergistic estimation of AGB using Landsat 8 OLI- and Sentinel-1A-derived products demonstrated better performance.

Overall, the AGB ground reference data derived from the Tent_ASO_BP model with integrated variables from Landsat 8 OLI and Sentinel-1A as inputs achieved satisfactory biomass estimation precision. The proposed modeling efforts are expected to greatly facilitate future large-scale forest growth monitoring and carbon stock assessment.