Inversion of Forest Aboveground Biomass in Regions with Complex Terrain Based on PolSAR Data and a Machine Learning Model: Radiometric Terrain Correction Assessment

Abstract

The accurate estimation of forest aboveground biomass (AGB) in areas with complex terrain is very important for quantifying the carbon sequestration capacity of forest ecosystems and studying the regional or global carbon cycle. In our previous research, we proposed the radiometric terrain correction (RTC) process for introducing normalized correction factors, which has strong effectiveness and robustness in terms of the backscattering coefficient of polarimetric synthetic aperture radar (PolSAR) data and the monadic model. However, the impact of RTC on the correctness of feature extraction and the performance of regression models requires further exploration in the retrieval of forest AGB based on a machine learning multiple regression model. In this study, based on PolSAR data provided by ALOS-2, 117 feature variables were accurately extracted using the RTC process, and then Boruta and recursive feature elimination with cross-validation (RFECV) algorithms were used to perform multi-step feature selection. Finally, 10 machine learning regression models and the Optuna algorithm were used to evaluate the effectiveness and robustness of RTC in improving the quality of the PolSAR feature set and the performance of the regression models. The results revealed that, compared with the situation without RTC treatment, RTC can effectively and robustly improve the accuracy of PolSAR features (the Pearson correlation R between the PolSAR features and measured forest AGB increased by 0.26 on average) and the performance of regression models (the coefficient of determination R² increased by 0.14 on average, and the rRMSE decreased by 4.20% on average), but there is a certain degree of overcorrection in the RTC process. In addition, in situations where the data exhibit linear relationships, linear models remain a powerful and practical choice due to their efficient and stable characteristics. For example, the optimal regression model in this study is the Bayesian Ridge linear regression model (R² = 0.82, rRMSE = 18.06%).

1. Introduction

Forest aboveground biomass (AGB) is a core parameter for monitoring forest ecosystems and assessing carbon sink capacity. Therefore, it plays an important role in quantifying the carbon sequestration capacity of forest ecosystems and studying regional and global carbon cycles [1,2]. However, it is challenging to accurately estimate forest AGB, especially in regions with complex and variable topography, which increases the uncertainty of assessments of forest ecosystem productivity. Accuracy is important for improving forest resource status surveys [3,4]. Remote sensing technology, including optical remote sensing, microwave remote sensing, and light detection and ranging (LiDAR), has become the most practical way to estimate forest AGB because it is long-term, efficient, and economical for acquiring large-scale spatial distribution information about forests [5,6,7]. Among the potential techniques, synthetic aperture radar (SAR) has been widely used for the dynamic monitoring of forest AGB in large-scale areas due to its good performance in all weathers, sustainable observation, and cloud penetration [8,9,10]. Although multi-source remote sensing data can significantly improve the inversion accuracy of forest AGB [11,12,13,14,15], it is still crucial to explore AGB estimation methods relying only on SAR data. Especially in areas where optical remote sensing is limited due to frequent cloud and fog, SAR becomes an irreplaceable and effective tool [16].

At present, there are four common technologies for forest AGB based on SAR data: interferometric SAR (InSAR) [17], polarimetric SAR (PolSAR) [18,19], polarimetric interferometric SAR (PolInSAR) [20,21,22], and tomographic SAR (TomoSAR) [23,24]. PolSAR technology is widely used in large-scale inversion of forest AGB because it does not rely on interference technology and does not need to meet the strict requirements of repeated orbit observation and accurate time synchronization. It has more relaxed data acquisition conditions and simplicity of data collection. In addition, the modeling methods for estimating forest AGB based on PolSAR data can be divided into scattering mechanism methods [25], machine learning methods [26,27,28], and deep learning methods [29,30,31]. Scattering mechanism methods (such as the water cloud model (WCM)) [32,33] can be useful because with a simplified physical model it is difficult to describe the real scattering characteristics of a forest with a complex structure. Therefore, the universal applicability of this kind of model may be poor in a forest area with strong heterogeneity. For deep learning methods, a large amount of labeled data are required for the model training. However, the actual measurement of forest biomass is often costly and difficult, and so such models are relatively poor in interpretation and prone to overfitting. Therefore, this study is based on PolSAR data, using a flexible and more adaptable machine learning method to estimate forest AGB.

The process of inversion of forest AGB by the machine learning method based on PolSAR data includes three key stages: firstly, accurate extraction and derivation of PolSAR feature parameters; secondly, selection of the optimal feature subset; and, finally, selection of the most suitable regression model and algorithm. Accurate extraction of PolSAR features allows us to mine the information contained in the original data and generate new features based on physical mechanisms or statistical methods. This enhances the ability of the statistical model to reflect the complex relationship between forest AGB and PolSAR features [34]. Backscattering coefficient and polarization decomposition parameters are the most commonly used PolSAR features in estimating forest AGB, with their extraction theories having been extensively studied in literature [35,36,37]. However, due to the side-view imaging principle of the SAR system and the characteristics of tilt distance being inferred based on echo delay, the actual topographic relief and other factors can cause serious interference to the scattered echo information received by SAR. This can limit the accuracy and reliability of the features extracted from PolSAR data [38]. Therefore, fully considering the influence of terrain and implementing effective terrain correction is essential to improving the accuracy of PolSAR features. At present, the radiometric terrain correction (RTC) process designed for PolSAR data is a relatively complete solution for terrain correction [39], which includes polarization orientation angle correction (POAC), effective scattering area correction (ESAC), and angular variation effect correction (AVEC). We introduce normalized ESAC factors on the basis of RTC to further improve the RTC process [40]. Although the RTC method has shown excellent performance in the correction of backscattering coefficient and its effectiveness and robustness have been widely recognized [41,42], there is a lack of in-depth research on the correction effect and quantitative analysis of the polarization decomposition parameters of PolSAR data by RTC. In addition, the effect of the fusion of the RTC process with different machine learning algorithms and potential improvements need to be discussed further.

Feature selection is essential to achieve robust and high-precision estimation of forest AGB based on PolSAR data [43,44,45]. The data and features determine the upper limit of machine learning, while models and algorithms only approach this upper limit [46]. The purpose of selecting the optimal feature subset is to remove invalid, redundant, and interference-producing variables from many potential feature variables, and solve the multicollinearity problem between the features synchronously. The aim is to select the feature set with the greatest information value and prediction efficiency so as to improve the efficiency and accuracy of the model construction [47,48,49,50]. At present, the methods of feature selection are mainly divided into three categories: filter, wrapper, and embedded algorithms [51,52]. In addition, a hybrid feature selection algorithm combining the advantages of multiple feature selection methods through multi-step screening has also been proposed [53,54,55].

The selection of appropriate statistical models can help to build models that accurately reflect the complex relationship between AGB and remote sensing data, thus improving the reliability and accuracy of carbon storage assessment. At present, according to the different construction methods of statistical models and the characteristics of parameter assumptions, regression models based on machine learning can be divided into two categories: parametric models and non-parametric models [5,56]. Non-parametric models can flexibly adapt to various complex data structures without pre-setting the specific form of the model, and thus show excellent adaptability and robustness. However, the parametric model is still a competitive regression analysis tool because of its clear explanations and high computational efficiency [57]. In addition, in the inversion of forest AGB based on PolSAR characteristics, a widely accepted view is that forest AGB has the best correlation with the backscattering coefficient of HV polarization channels, but there is no direct and simple linear relationship between the two; in fact, a logarithmic function may be the more appropriate relationship [58]. Moreover, our previous research focused on the potential of RTC processes to improve the accuracy of forest AGB inversion based on univariate models [40]. Therefore, this study will not focus on the construction of monadic linear or nonlinear models, but instead will seek multivariate and high-level nonlinear model structures that can capture this complex relationship.

This research aims to achieve three key tasks: (1) quantitatively evaluate the ability of the RTC process to enhance the correlation between forest AGB and PolSAR features (backscattering coefficient, polarization decomposition parameters, and related derived features); (2) explore the effect of RTC on polarization decomposition parameters and the underlying mechanism of action; and (3) evaluate the potential of the RTC process to improve the performance of forest AGB inversion based on machine learning models, and identify the optimal model by comparing the performance of different models.

2. Materials

2.1. Study Area

The study area is located in Saihanba Forest Farm (116°53′~117°43′E, 41°55′~42°35′N, Figure 1) in the border area between the northern mountains of Weichang County, Hebei Province, China, with an altitude of 1010 to 1940 m. The study area covers a wide range of vegetation types, including deciduous coniferous forest, evergreen coniferous forest, broad-leaved forest, scrubland, and grassland. Moreover, the topography of the study area is complex, with the forest cover being a combination of flat woodlands and steep hillsides with slopes up to 55.47°. This variety of terrain conditions provides an ideal natural laboratory for conducting comparative studies of complex terrain. In addition, the study area has a large-scale artificial forest, covering an area of about 93,461 hectares, and its large area of artificial forest is convenient for scientific researchers to collect real and rich test data on the spot.

2.2. PolSAR Data and Pre-Processing

The PolSAR data used in this study are L-band full-polarization single-look complex (SLC) data (

Table 1. The main parameters of the PolSAR data.

Parameter	Value	Parameter	Value
Data formats and processing level	CEOS level 1.1	Observation mode	HBQ
Observation date of scene center	27.8054°	Radar wavelength	0.2424525 m
Length of range direction	49.8 km	Range resolution	2.860844 m
Length of azimuth direction	69.3 km	Azimuth resolution	2.642742 m
Orbit direction	Ascending	Range pixel	8392
Observation direction	Right looking	Azimuth pixel	26,105

Note: HBQ is High-sensitive mode Full (Quad.) polarimetry.

) obtained by the Advanced Land Observing Satellite-2 (ALOS-2) (developed by the Japan Aerospace Exploration Agency (JAXA)) on 11 July, 25 July, and 8 August 2020. As the only platform currently equipped with a satellite-borne L-band SAR sensor, ALOS-2’s four-polarization observation mode can significantly enhance its ability to capture complex forest structures and scattering mechanisms. Moreover, compared with C-band, L-band has stronger penetration and can more accurately reveal the details of the forest internal structure, which is conducive to improving the accuracy of forest AGB.

In this study, the Sentinel Application Platform (SNAP_v9.0.0) software developed by the European Space Agency completed the pre-processing steps for the PolSAR data: cross-channel SNR correction, Faraday rotation correction, and radiometric calibration for complex data (Equation (2), derived from radiometric calibration of real data; Equation (1)) [35,40], multi-look averaging of 4 × 9 (range × azimuth), refined LEE filter, and range doppler terrain correction). The ensemble average size for the covariance matrix is crucial, especially when studying natural vegetation. For natural vegetation, when the reflection symmetry condition remains unchanged, the volume scattering power induced by vegetation will improve with the increase in the average ensemble size [59]. The purpose of choosing 4 × 9 (range × azimuth) as the core parameter of multi-look processing is to process the pixels of the PolSAR data into square pixels as close as possible to the reference DEM data and the measured plot.

$${\sigma }_{slc}^0=10\cdot {\log }_{10}〈I^2+Q^2〉+C{F}_1-A,$$

(1)

where ${\sigma }_{slc}^0$ is the backscattering coefficient after radiation calibration; I and Q represent the real and imaginary parts of SLC 1.1 products, respectively; CF₁ is −83.0 dB; and A is 32.0 dB.

$$I_{cal}=\frac{I}{{10}^{5.75}};Q_{cal}=\frac{Q}{{10}^{5.75}},$$

(2)

where I_cal and Q_cal are the real and imaginary parts of SLC 1.1 products after radiation calibration.

2.3. Ground-Measured Forest AGB

The measured forest AGB data were collected through a field survey conducted at Saihanba Forest Farm in August 2020 for the training and validation of the forest AGB inversion model. Through screening the key stand characteristics of the plots, including tree species composition, understory environment, topographic slope, and stand thinning density, 132 rhomboid plots (Figure 1) with an area of 0.06 ha (24.49 m × 24.49 m) were established and investigated, including 28 plots distributed on a 1.0 km grid and 104 plots with a relatively uniform distribution. All plots were located using differential GPS (Unistrong RTK-G10, Unistrong, Bei**g, China), and plot information such as plot type, plant number, and slope, as well as individual information of trees with a diameter at breast height (DBH) ≥ 5.0 cm, was recorded, including DBH, tree height, and tree species. Finally, the biomass per plant was calculated separately according to the allometric equation of different tree species (

Table 2. Allometric growth equation of different tree species.

Tree Species	Allometric Equation	Ref.
Larix principis-rupprechtii	W = 0.1431 · DBH2.2193	[60]
Betula platyphylla	W = 0.0330 · DBH2.9314	[61]
Pinus sylvestris var. mongolica	W = 0.0930 · DBH2.3429	[62]
Pinus tabuliformis	W = 0.0520 · DBH2.5830	[63]
Acer truncatum	W = 0.1260 · DBH2.3830	[64]

Note: W is the AGB (kg) per plant, and DBH is the diameter at breast height (cm).

); then the total forest AGB of the plot was obtained by summing (Equation (3)) [40] all the biomass per plant.

Table 3. Sample plot of AGB statistics.

Count	Min (t/ha)	Max (t/ha)	Mean (t/ha)	SD (t/ha)
132	36.39	166.58	103.88	32.31

Note: SD is standard deviation; Count is the number of sample plots.

provides statistical details of the field data for this study.

$$AG{B}_j={\displaystyle {\sum }_{i=1}^m{W}_i}/1000\times S,$$

(3)

where AGB_j is the AGB (t/ha) of the jth sample plot, m is the number of trees in the plot, W_i is the AGB of the ith tree (kg), and S is the area of the sample plot (S = 0.06 ha).

2.4. SRTM DEM

Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM) data with a resolution of 30 m were used to complement the PolSAR data to implement geocoding and calculate RTC-related factors.

3. Methods

The processing framework of this study (Figure 2) includes the following steps: (1) radiometric terrain correction (RTC) of PolSAR data, including polarization orientation angle correction (POAC), effective scattering area correction (ESAC), and angular variation effect correction (AVEC); (2) feature extraction and feature derivation based on PolSAR data, including backscattering coefficient, polarization decomposition parameters, and their related derived feature variables; and (3) feature selection and model training. All machine learning algorithms were executed in the Python 3.10 development environment. The programs that involved partitioning the dataset used the default parameter (test_size = 0.25), which was 75% for the training set and 25% for the test set. The following sections present the necessary theories and related formulas for each stage in the flowchart in more detail.

3.1. Radiometric Terrain Correction for PolSAR

3.1.1. Polarization Orientation Angle Correction

When there is topographic relief in the target area, the polarization orientation angle (POA)—that is, the angle between the long axis of the polarization ellipse and the horizontal axis of the ground plane—will be biased due to the slope (mainly the azimuth slope). In order to compensate for this deviation of POA caused by terrain, the POA extracted by the circular polarization method [39,65] is used first, and then the three-dimensional polarization covariance (C3) matrix is compensated (Equation (4)).

$$\begin{array}{l}{C}_{3\_POAC}=U_{3\left(\eta \right)}{C}_3{U}_{3\left(\eta \right)}^{-1}\\ U_{3\left(\eta \right)}=\frac{1}{2}\left[\begin{array}{ccc}1+\cos \left(2\eta \right)& \sqrt{2}\sin \left(2\eta \right)& 1-\cos \left(2\eta \right)\\ -\sqrt{2}\sin \left(2\eta \right)& 2\cos \left(2\eta \right)& \sqrt{2}\sin \left(2\eta \right)\\ 1-\cos \left(2\eta \right)& -\sqrt{2}\sin \left(2\eta \right)& 1+\cos \left(2\eta \right)\end{array}\right]\\ \eta =\frac{1}{4}\left[-1\times Arg\left(〈\frac{-4\mathrm{Re}〈\left(S_{HH}-S_{VV}\right)S_{HV}^{*}〉}{-〈{\left|S_{HH}-S_{VV}\right|}^2〉+4〈{\left|S_{HV}\right|}^2〉}〉\right)+\pi \right]\end{array},$$

(4)

where C₃ and C_{3_POAC} are the C₃ matrix before and after polarization orientation angle correction (POAC), respectively; η is the polarization orientation angle (POA); Arg(·) is a phase function; and when η > π/4, η = η − π/2.

3.1.2. Effective Scattering Area Correction

The effective scattering area is the projected area of the actual ground element on the isophase surface. However, due to the difference in terrain slope, ground units of the same SAR data in the geographical coordinate system correspond to different effective scattering areas. Therefore, effective scattering area correction (ESAC) takes flat land as a reference to correct the area effect of SAR data under different slope conditions so as to eliminate the influence of terrain slope on scattering information [39,40]. The C₃ matrix can be corrected using Equation (5):

$$C_{3\_ESAC}=C_3\cdot \frac{\sin \left({\theta }_{loc}\right)}{\sin \left({\theta }_{ref}\right)},$$

(5)

where C_{3_ESAC} is the C₃ matrix after ESAC; C₃ refers to the C₃ matrix without ESAC; θ_loc is the local incidence angle; and θ_ref is the radar incidence angle.

3.1.3. Angular Variation Effect Correction

The angular variation effect (AVE) refers to the phenomenon whereby the local incidence angle changes due to a change in terrain, which leads to a change in the scattering mechanism of the target object (such as forest). The angular variation effect correction (AVEC) is typically based on the cosine model (Equation (6)) to correct forest-covered areas. In this study, the forest cover area was first extracted by Freeman3 decomposition and Wishart unsupervised classification [66], and then the optimal n value matrix of the C₃ matrix (different elements correspond to different optimal n values) was obtained by calculating the optimal n value corresponding to the main diagonal elements of the C₃ matrix. For each element of the three principal diagonal elements, the optimal n value was first determined by iterating the minimum value of the absolute value of f(n) when n takes different values [39].

$$\begin{array}{c}{C}_{3\_AVEC}=C_3\cdot {\left(\frac{\cos {\theta }_{ref}}{\cos {\theta }_{loc}}\right)}^{n_{pq}}=C_3\odot \left[\begin{array}{ccc}k\left(n_{HH}\right)& k\left(\frac{n_{HH}+n_{HV}}{2}\right)& k\left(\frac{n_{HH}+n_{VV}}{2}\right)\\ k\left(\frac{n_{HH}+n_{HV}}{2}\right)& k\left(n_{HV}\right)& k\left(\frac{n_{HV}+n_{VV}}{2}\right)\\ k\left(\frac{n_{HH}+n_{VV}}{2}\right)& k\left(\frac{n_{HV}+n_{VV}}{2}\right)& k\left(n_{VV}\right)\end{array}\right]\\ n_{pq}=argmin\left\{f_{\left(n\right)}\right\}=argmin\left\{\left|\rho \left({\theta }_{loc},C_{3\_AVEC\_pq}\right)\right|\right\}\end{array},$$

(6)

where (p, q) = (H, V) is the polarization mode of incidence and scattering of electromagnetic waves; n_pq is the optimal n value of different polarization channels; ⊙ is the Hadamard product; argmin{·} is the operation that takes the minimum value; f(n) is the absolute value of the correlation coefficient between the local incidence angle θ_loc and C_{3_AVEC_pq}; C_{3_AVEC_pq} is the decibel state of the intensity value after AVEC is applied to one of the principal diagonal elements of the C₃ matrix at a certain value of n; and ρ(·) is the correlation coefficient function.

3.2. Feature Extraction and Feature Derivation of PolSAR

All 118 PolSAR features extracted in this study were based on the Polarimetric SAR Data Processing and Educational Toolbox (PolSARpro_v5.1.1) and the Sentinel Application Platform (SNAP_v9.0.0).

3.2.1. Backscattering Coefficient Features of PolSAR

This study extracted the intensity values (${\sigma }_{HH}^0$, ${\sigma }_{HV}^0$, and ${\sigma }_{VV}^0$) and decibel state values (${\sigma }_{HH}^0$_db, ${\sigma }_{HV}^0$_db, and ${\sigma }_{VV}^0$_db) of the three main diagonal elements of the C₃ matrix as the six original backscatter coefficient features, and then calculated 20 classical derived features (

Table 4. The features derived from the backscattering coefficients extracted from the PolSAR data.

Features	Symbol	Equation	Symbol (dB)	Source
Span	Span	${\sigma }_{HH}^0+2{\sigma }_{HV}^0+{\sigma }_{VV}^0$	Span_db	[67]
Co-Pol HH/VV Ratio	R_HH/VV	${\sigma }_{HH}^0/{\sigma }_{VV}^0$	R_HH/VV_db	[67]
Cross-Pol HH/HV Ratio	R_HH/HV	${\sigma }_{HH}^0/{\sigma }_{HV}^0$	R_HH/HV_db	[67]
Cross-Pol VV/VH Ratio	R_VV/VH	${\sigma }_{VV}^0/{\sigma }_{VH}^0$	R_VV/VH_db	[68]
Biomass Index	BMI	$\left({\sigma }_{HH}^0+{\sigma }_{VV}^0\right)/2$	BMI_db	[69]
Volume Scattering Index	VSI	${\sigma }_{HV}^0/\left({\sigma }_{HV}^0+BMI\right)$	VSI_db	[69]
Canopy Structure Index	CSI	${\sigma }_{VV}^0/\left({\sigma }_{HH}^0+{\sigma }_{VV}^0\right)$	CSI_db	[69]
Radar Vegetation Index	RVI	$8{\sigma }_{HV}^0/\left({\sigma }_{HH}^0+2{\sigma }_{HV}^0+{\sigma }_{VV}^0\right)$	RVI_db	[70]
Radar Forest Degradation Index	RFDI	$\left({\sigma }_{HH}^0-{\sigma }_{HV}^0\right)/\left({\sigma }_{HH}^0+{\sigma }_{HV}^0\right)$	RFDI_db	[71]
modified RFDI	mRFDI	$\left({\sigma }_{VV}^0-{\sigma }_{HV}^0\right)/\left({\sigma }_{VV}^0+{\sigma }_{HV}^0\right)$	mRFDI_db	[72]

Note: The backscattering coefficient ${\sigma }_{pq}^0$ in the table is the backscattering intensity in linear units.

).

3.2.2. Polarization Decomposition Features of PolSAR

In this study, 72 original polarization decomposition-related features in linear units and decibel units were extracted based on 12 polarization decomposition methods (

Table 5. The polarization decomposition features extracted from the PolSAR data.

Decomposition Methods	Abbreviation	Symbol	Source
Freeman two-component	FRE2	F2V, F2O	[73]
Freeman three-component	FRE3	F3V, F3D, F3O	[74]
Yamaguchi three-component	YAM3	Y3V, Y3D, Y3O	[75]
Yamaguchi four-component	YAM4	Y4V, Y4D, Y4O, Y4H	[76]
Van Zyl three-component	VAZ3	V3V, V3D, V3O	[77]
An and Yang three-component	ANY3	A3V, A3D, A3O	[78]
Model-free three-component	MF3CF	M3V, M3D, M3O	[79]
Model-free four-component	MF4CF	M4V, M4D, M4O, M4H	[80]
Cloude three-component	CLD3	C3V, C3D, C3O	[81]
Krogager three-component	KRO3	K3S, K3D, K3H	[82]
H-A-alpha	HAα	H, A, α	[83]
Pauli three-component	PAU3	P3D, P3O	[84]

Note: In this table, the symbol of each component of the different polarization decomposition methods (except H-A-alpha) has a consistent marking method: the first character (letter) is the first letter of the abbreviation of the polarization decomposition method, the second character (number) is the number of its components, and the third character (letter) is the first letter of the abbreviation of the different component types, including Volume scattering component (Vol), Surface scattering component (Odd), Double-bounce scattering component (Dbl), Helix scattering component (Hel), and Sphere scattering component (Sph). In addition, the decibel symbol of each polarization decomposition component adds a suffix (“_db”) to the symbol of its intensity unit, which is not listed in this table. For example, the decibel symbol of F3V is F3V_db.

), and then we calculated 20 derived features (

Table 6. The features derived from the polarization decomposition extracted from the PolSAR data.

Symbol	Equation	Symbol (dB)	Source
F2R_1	$F2V/\left(F2O\right)$	F2R_1_db	[85]
F3R_1	$F3V/\left(F3D+F3O\right)$	F3R_1_db	[85]
A3R_1	$A3V/\left(A3D+A3O\right)$	A3R_1_db	[85]
V3R_1	$V3V/\left(V3D+V3O\right)$	V3R_1_db	[85]
Y3R_1	$Y3V/\left(Y3D+Y3O\right)$	Y3R_1_db	[85]
F2R_2	$F2V/\left(F2V+F2O\right)$	F2R_2_db	[86]
F3R_2	$F3V/\left(F3V+F3D+F3O\right)$	F3R_2_db	[86]
A3R_2	$A3V/\left(A3V+A3D+A3O\right)$	A3R_2_db	[86]
V3R_2	$V3V/\left(V3V+V3D+V3O\right)$	V3R_2_db	[86]
Y3R_2	$Y3V/\left(Y3V+Y3D+Y3O\right)$	Y3R_2_db	[86]

Note: The meanings of the symbols in the equations in this table are shown in Table 5.

). The Volume scattering component (Vol) of the Pauli three-component (PAU3) in this study is exactly the same as σ⁰_HV, so only the Surface scattering component (Odd) and the Double-bounce scattering component (Dbl) of the PAU3 decomposition method were used.

3.3. Forest AGB Regression Modeling Algorithms and Model Evaluation

In this study, 10 multiple regression models (

Table 7. The 10 regression models used in this study.

Model	Python Package	Module	Estimator	Category
Ridge	scikit-learn (v1.4.2)	sklearn.linear_model	Ridge	linear
Lasso	scikit-learn (v1.4.2)	sklearn.linear_model	Lasso	linear
ElasticNet	scikit-learn (v1.4.2)	sklearn.linear_model	ElasticNet	linear
BysRidge	scikit-learn (v1.4.2)	sklearn.linear_model	BayesianRidge	linear
ARD	scikit-learn (v1.4.2)	sklearn.linear_model	ARDRegression	linear
RF	scikit-learn (v1.4.2)	sklearn.ensemble	RandomForestRegressor	non-parametric
ET	scikit-learn (v1.4.2)	sklearn.ensemble	ExtraTreesRegressor	non-parametric
AdaBoost	scikit-learn (v1.4.2)	sklearn.ensemble	AdaBoostRegressor	non-parametric
XGBoost	xgboost (v2.0.3)	xgboost	XGBRegressor	non-parametric
CatBoost	catboost (v1.2.3)	catboost	CatBoostRegressor	non-parametric

) were selected to evaluate the effectiveness and robustness of RTC on the performance of different machine learning regression models. Multiple regression models are statistical and machine learning methods used to analyze the relationship between a continuous dependent variable (response variable) and two or more independent variables (feature variables), which can be divided into multiple linear regression models and non-parametric regression models. Non-parametric regression models do not depend on the form of the model, while the multiple linear regression model has an explicit parametric form (Equation (7)) [67]:

$$Y={\alpha }_0+{\alpha }_1{X}_1+{\alpha }_2{X}_2+\cdots +{\alpha }_n{X}_n+\epsilon ,$$

(7)

where Y is the target variable of the prediction (decibels of forest AGB); X₁, X₂, …, X_n are the predictor variables; α₀ is a constant; α₁, α₂, …, α_2n are the regression coefficients associated with the corresponding variables; n is the number of the predictor variables; and ε is the error term.

The machine learning process in this study included the following steps (implementation details shown in

Table 8. The Python packages on which the implementation of two-step feature selection depends.

Algorithm	Boruta	RFECV	RF	Optuna
package	Boruta (v0.3)	scikit-learn (v1.4.2)	scikit-learn (v1.4.2)	Optuna (v3.6.0)

): (1) The random forest regression model (RF) based on the default parameters served as the base model for Boruta’s algorithm [87] and performed the initial feature filtering on the entire dataset instead of the training set. The goal is to identify and remove features that are irrelevant or weakly related to the target variable. (2) First, the Optuna algorithm (n_trials = 500, including 10-fold cross-validation) was used to obtain the optimal hyperparameters of 10 regression models on the feature set after the initial feature selection. Then, the 10 regression models based on the optimal hyperparameters corresponding to different models were used as the base model for recursive feature elimination with cross-validation (RFECV) feature selection, and the feature selection of the second step was completed. The aim is to iteratively search the optimal feature subset of the base model. (3) Again, the Optuna algorithm was used to obtain and output the corresponding evaluation indicators (based on the test set, i.e., 25% of all sample plots) of optimal performance of 10 regression models on the feature set after two-step feature selection. These evaluation indicators included the coefficient of determination (R², Equation (8)), the root mean square error (RMSE, Equation (9)), and the relative root mean square error (rRMSE, Equation (10)) [18,27,42]. This was repeated 10 times and we calculated the average of the above indicators. This is because the performance of the model largely depends on the selected hyperparameters; however, the values of these hyperparameters are often sensitive to the dataset, so the hyperparameters of Optuna need to be tuned separately for different feature sets. In addition, since random states cannot be completely controlled, they were not fixed in this study. This means that the initial state of model training and possible random factors in each experiment followed a natural (random) process, which was more in line with the real situation. However, through cross-validation and repeated tests, randomness can be reduced and a more stable and representative evaluation index can be obtained.

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

(8)

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

(9)

$$rRMSE=\left(RMSE/\overline{y}\right)\times 100\%$$

(10)

where N is the sample size, $y_i$ represents the true observed value, $\overline{y}$ represents the average of the true observed value, and ${\widehat{y}}_i$ represents the predicted value.

4. Results

4.1. The Impact of Radiometric Terrain Correction on PolSAR Features

In order to quantitatively analyze the impact of the radiometric terrain correction (RTC) process on the correlation between forest AGB and PolSAR features, we first set the non-RTC (NRTC) data as the blank control and the RTC data as the test data. The core correction factors of the RTC process are shown in our previous study [40]. Secondly, we extracted a total of 117 effective PolSAR features (as mRFDI_db is invalid) based on NRTC and RTC data. Subsequently, Pearson correlation coefficients (R) were calculated separately between forest AGB and each of the PolSAR features (

Table A1. Pearson correlation coefficient (R) between forest AGB and the 118 PolSAR features based on the data with radiometric terrain correction (RTC) and non-RTC data (NRTC), obtained from 11 July, 25 July, and 8 August 2020.

ID	Features	NRTC/RTC			Significance
		R_200711	R_200725	R_200808
1	${\mathsf{\sigma }}_{\mathrm{HH}}^0$	0.12/0.58	0.16/0.63	0.14/0.60	-/**
2	${\mathsf{\sigma }}_{\mathrm{HV}}^0$	0.25/0.70	0.28/0.73	0.27/0.71	*/**
3	${\mathsf{\sigma }}_{\mathrm{VV}}^0$	0.15/0.61	0.18/0.66	0.17/0.62	-/**
4	Span	0.18/0.66	0.21/0.71	0.20/0.69	-/**
5	R_HH/VV	−0.09/−0.09	−0.10/−0.10	−0.10/−0.10	-/-
6	R_HH/HV	−0.74/−0.73	−0.77/−0.76	−0.76/−0.75	**/**
7	R_VV/VH	−0.70/−0.69	−0.74/−0.73	−0.72/−0.70	**/**
8	BMI	0.14/0.62	0.17/0.66	0.15/0.65	-/**
9	VSI	0.70/0.67	0.73/0.71	0.71/0.68	**/**
10	RVI	0.70/0.67	0.73/0.71	0.71/0.68	**/**
11	CSI	0.091/0.089	0.096/0.094	0.096/0.095	-/-
12	RFDI	−0.65/−0.64	−0.68/−0.66	−0.67/−0.65	**/**
13	mRFDI	−0.63/−0.61	−0.66/−0.64	−0.65/−0.62	**/**
14	F2O	0.14/0.51	0.19/0.56	0.16/0.53	-/**
15	F2V	0.14/0.51	0.17/0.54	0.14/0.52	-/**
16	F3D	0.15/0.30	0.17/0.31	0.19/0.34	-/*
17	F3O	0.04/0.11	0.05/0.14	0.03/0.10	-/-
18	F3V	0.30/0.72	0.32/0.74	0.31/0.74	*/**
19	Y3D	0.22/0.34	0.26/0.36	0.24/0.35	-/*
20	Y3O	−0.03/0.10	−0.04/0.15	−0.04/0.11	-/-
21	Y3V	0.25/0.69	0.29/0.73	0.28/0.70	-/**
22	Y4D	−0.04/0.21	−0.05/0.23	−0.04/0.21	-/-
23	Y4O	−0.24/−0.44	−0.25/−0.44	−0.27/−0.45	-/*
24	Y4V	0.31/0.71	0.30/0.71	0.30/0.70	*/**
25	Y4H	0.08/0.15	0.09/0.16	0.09/0.18	-/-
26	M3D	0.21/0.49	0.23/0.51	0.20/0.51	-/**
27	M3O	0.01/0.25	0.01/0.29	0.01/0.28	-/-
28	M3V	0.32/0.72	0.30/0.71	0.35/0.74	*/**
29	M4D	0.25/0.60	0.24/0.58	0.22/0.57	-/**
30	M4O	0.01/0.28	0.01/0.28	0.01/0.30	-/-
31	M4V	0.29/0.71	0.30/0.72	0.32/0.74	*/**
32	M4H	0.06/0.13	0.06/0.12	0.09/0.15	-/-
33	A3D	0.22/0.45	0.20/0.43	0.18/0.40	-/*
34	A3O	0.03/0.32	0.03/0.31	0.03/0.30	-/-
35	A3V	0.31/0.64	0.29/0.63	0.30/0.63	-/**
36	V3D	0.20/0.46	0.22/0.49	0.22/0.51	-/*
37	V3O	0.15/0.54	0.16/0.54	0.18/0.55	-/**
38	V3V	0.32/0.68	0.31/0.66	0.31/0.65	*/**
39	H	0.61/0.62	0.61/0.63	0.60/0.63	**/**
40	A	−0.07/−0.23	−0.06/−0.23	−0.07/−0.25	-/-
41	α	0.55/0.51	0.57/0.55	0.57/0.56	**/**
42	C3D	0.24/0.48	0.24/0.49	0.26/0.50	-/*
43	C3O	0.30/0.60	0.32/0.63	0.31/0.60	-/**
44	C3V	−0.01/0.15	−0.02/0.19	−0.02/0.20	-/-
45	K3S	0.18/0.52	0.20/0.55	0.21/0.55	-/**
46	K3D	0.31/0.70	0.30/0.68	0.30/0.67	*/**
47	K3H	0.05/0.08	0.05/0.08	0.05/0.07	-/-
48	P3D	0.20/0.60	0.22/0.63	0.22/0.65	-/**
49	P3O	0.13/0.52	0.14/0.54	0.15/0.55	-/**
50	A3R_1	0.60/0.55	0.60/0.56	0.61/0.58	**/**
51	A3R_2	0.62/0.59	0.62/0.58	0.63/0.60	**/**
52	F2R_1	−0.31/−0.12	−0.33/−0.13	−0.37/−0.15	-/-
53	F2R_2	−0.32/−0.15	−0.34/−0.15	−0.38/−0.17	-/-
54	F3R_1	0.61/0.58	0.63/0.59	0.59/0.55	**/**
55	F3R_2	0.69/0.60	0.70/0.61	0.71/0.61	**/**
56	V3R_1	0.70/0.66	0.71/0.68	0.69/0.66	**/**
57	V3R_2	0.72/0.70	0.73/0.70	0.71/0.69	**/**
58	Y3R_1	0.65/0.61	0.65/0.62	0.64/0.60	**/**
59	Y3R_2	0.70/0.64	0.72/0.64	0.70/0.63	**/**
60	${\mathsf{\sigma }}_{\mathrm{HH}}^0\text{\_}\mathrm{db}$	0.30/0.70	0.31/0.72	0.31/0.71	*/**
61	${\mathsf{\sigma }}_{\mathrm{HV}}^0\text{\_}\mathrm{db}$	0.50/0.81	0.54/0.86	0.52/0.82	**/**
62	${\mathsf{\sigma }}_{\mathrm{VV}}^0\text{\_}\mathrm{db}$	0.32/0.72	0.34/0.75	0.34/0.73	*/**
63	Span_db	0.38/0.79	0.40/0.82	0.40/0.80	*/**
64	R_HH/VV_db	−0.01/−0.01	−0.098/−0.096	−0.10/−0.10	-/-
65	R_HH/HV_db	−0.57/−0.65	−0.71/−0.69	−0.70/−0.67	**/**
66	R_VV/VH_db	−0.62/−0.60	−0.66/−0.65	−0.64/−0.63	**/**
67	BMI_db	0.30/0.73	0.32/0.76	0.30/0.75	*/**
68	VSI_db	0.75/0.74	0.77/0.75	0.78/0.76	**/**
69	RVI_db	0.75/0.74	0.77/0.75	0.78/0.76	**/**
70	CSI_db	0.10/0.09	0.099/0.097	0.10/0.10	-/-
71	RFDI_db	−0.62/−0.60	−0.64/−0.61	−0.65/−0.62	**/**
72	mRFDI_db	NaN/NaN	NaN/NaN	NaN/NaN	-/-
73	F2O_db	0.36/0.66	0.40/0.69	0.38/0.67	*/**
74	F2V_db	0.23/0.53	0.25/0.56	0.25/0.54	-/**
75	F3D_db	0.35/0.45	0.37/0.46	0.38/0.48	*/*
76	F3O_db	−0.05/0.14	−0.06/0.16	−0.04/0.11	-/-
77	F3V_db	0.45/0.81	0.47/0.83	0.46/0.83	*/**
78	Y3D_db	0.36/0.47	0.39/0.49	0.37/0.47	*/*
79	Y3O_db	−0.03/0.10	−0.05/0.14	−0.05/0.12	-/-
80	Y3V_db	0.42/0.72	0.45/0.75	0.43/0.72	*/**
81	Y4D_db	−0.05/0.27	−0.07/0.29	−0.05/0.27	-/-
82	Y4O_db	−0.39/−0.65	−0.40/−0.65	−0.42/−0.67	*/**
83	Y4V_db	0.54/0.84	0.53/0.84	0.53/0.83	**/**
84	Y4H_db	0.16/0.17	0.18/0.19	0.18/0.20	-/-
85	M3D_db	0.39/0.60	0.40/0.62	0.38/0.62	*/**
86	M3O_db	0.04/0.32	0.05/0.35	0.05/0.34	-/*
87	M3V_db	0.55/0.83	0.54/0.82	0.58/0.85	**/**
88	M4D_db	0.41/0.63	0.40/0.62	0.39/0.60	*/**
89	M4O_db	0.05/0.35	0.05/0.35	0.05/0.37	-/*
90	M4V_db	0.55/0.82	0.55/0.83	0.56/0.85	**/**
91	M4H_db	0.12/0.17	0.12/0.16	0.15/0.18	-/-
92	A3D_db	0.39/0.57	0.38/0.56	0.35/0.52	*/**
93	A3O_db	0.05/0.38	0.05/0.37	0.04/0.35	-/*
94	A3V_db	0.50/0.79	0.49/0.79	0.49/0.80	**/**
95	V3D_db	0.38/0.55	0.41/0.57	0.41/0.59	*/**
96	V3O_db	0.28/0.63	0.29/0.63	0.30/0.64	-/**
97	V3V_db	0.46/0.79	0.46/0.78	0.45/0.76	**/**
98	H_db	0.62/0.63	0.62/0.64	0.61/0.64	**/**
99	A_db	−0.09/−0.06	0.08/−0.06	0.09/−0.08	-/-
100	α_db	0.54/0.50	0.56/0.53	0.56/0.54	**/**
101	C3D_db	0.33/0.57	0.33/0.58	0.35/0.59	*/**
102	C3O_db	0.37/0.70	0.40/0.72	0.38/0.70	*/**
103	C3V_db	0.02/0.23	0.04/0.26	0.04/0.27	-/-
104	K3S_db	0.31/0.58	0.34/0.61	0.35/0.61	*/**
105	K3D_db	0.47/0.78	0.46/0.76	0.46/0.75	*/**
106	K3H_db	0.09/0.09	0.089/0.092	0.09/0.09	-/-
107	P3D_db	0.35/0.63	0.38/0.65	0.38/0.67	*/**
108	P3O_db	0.27/0.55	0.28/0.57	0.28/0.58	-/**
109	A3R_1_db	0.65/0.59	0.65/0.60	0.66/0.63	**/**
110	A3R_2_db	0.66/0.63	0.66/0.62	0.67/0.64	**/**
111	F2R_1_db	−0.32/−0.15	−0.36/−0.16	−0.39/−0.18	*/-
112	F2R_2_db	−0.39/−0.19	−0.39/−0.19	−0.39/−0.19	*/-
113	F3R_1_db	0.64/0.62	0.66/0.63	0.61/0.60	**/**
114	F3R_2_db	0.71/0.64	0.72/0.65	0.73/0.65	**/**
115	V3R_1_db	0.73/0.69	0.75/0.71	0.73/0.70	**/**
116	V3R_2_db	0.75/0.72	0.76/0.72	0.75/0.71	**/**
117	Y3R_1_db	0.69/0.64	0.69/0.65	0.67/0.62	**/**
118	Y3R_2_db	0.73/0.67	0.75/0.67	0.73/0.66	**/**

Note: ** means extremely significant (p ≤ 0.01), * means significant (p ≤ 0.05), - means did not pass the significance test (p > 0.05). NaN means Not a Number, because mRFDI has a partial negative value, so it is impossible to calculate the feature of its decibel unit (mRFDI_db).

, Appendix B). Finally, in order to compare and visually present the changes in the R-values of the NRTC data and RTC data and, at the same time, reveal the relationship between various PolSAR features and forest AGB correlation strength after RTC, we sorted and plotted them according to the absolute R-values of the data (25 July 2020) after RTC (Figure 3).

The R-values of the RTC data increased to varying degrees in all (Figure 3a,b) but three original features (α, α_db, and A_db). The R-values of the 75 original features (with increased R-value) increased by an average of 0.26, and the largest R-value increase was 0.48 (${\sigma }_{HV}^0$), indicating that RTC processing can effectively improve the accuracy of most PolSAR features and the sensitivity of AGB inversion. However, among the original features, the three (H, H_db, and A) that belong to the same H-A-alpha (HAα) decomposition method as α, α_db, and A_db were relatively less affected by RTC. This shows that HAα decomposition can alleviate terrain disturbance. This is because the three parameters of the HAα decomposition are defined by a function of the ratio class of the three eigenvalues (the relative scattering power of different scattering mechanisms). Moreover, the R-values of the RTC data in all but four derived features (Span, Span_db, BMI, and BMI_db) (Figure 3c) were reduced to varying degrees. This suggests that the ratio operation (feature derived algorithm) of features can eliminate some multiplicative interference (including terrain interference). It also indicates that there is an overcorrection problem in the RTC process. Although the overcorrection problem of RTC is small, it causes an error of 0.05 in the mean value of the absolute value of the Pearson correlation coefficient between forest AGB and these 35 derived features.

In addition, the feature R_HH/HV (R = 0.76) has the largest R-value among all the features derived through ratio operations, but it is still 0.1 lower than ${\sigma }_{HV}^0$_db (R = 0.86), which has the highest R-value of all 117 PolSAR features. This suggests that, based only on the backscattering coefficient, the HV polarization mode carries more polarization information related to forest AGB. Although the ratio algorithm of features can mitigate the topographic effects, it will also weaken some AGB-related information. The second-largest R-value of all 117 PolSAR features is Y4V_db (R = 0.84), which indicates that the YAM4 decomposition method can decompose the scattering characteristics of the target into different scattering components relatively more accurately. Moreover, among the different components of the same polarization decomposition method, the Volume scattering component (Vol) usually has the largest R-value, indicating that the relationship between Vol and the scattering information in vegetation is relatively close. Furthermore, regardless of the decomposition method, the R-value of volume scattering in decibels is higher than that in linear units, which is consistent with the trend of previous studies [58].

4.2. The Impact of Radiometric Terrain Correction on Polarization Decomposition Component

In order to analyze the influence of different stages of RTC on the removal of topographic factors in the polarization decomposition component, we took the PolSAR data from 25 July 2020 as an example and plotted a scatter density plot between the decibel values of the three components of the Freeman three-decomposition method and the local incidence angle. As can be seen from Figure 4, there is a relatively obvious linear relationship between each component and the local incidence angle in the NRTC and POAC stages; but with the implementation of ESAC and AVEC, this linear relationship becomes significantly weaker, which is consistent with the results of our previous research regarding the backscattering coefficient [40]. The results show that RTC is also effective at removing the effect of terrain on the polarization decomposition.

In order to further analyze the specific mechanism of influence of different stages of RTC (POAC, ESAC, and AVEC) on each component of polarization decomposition, taking PolSAR data from 25 July 2020 as an example, we plotted the scatter density (Figure 5) of each component of Freeman3 decomposition at each stage of different RTC relative to the previous stage.

In the comparison stage between POAC and NRTC (Figure 5a,e,i), POAC caused the intensity values of almost all pixels of the Vol component (Figure 5a) to decrease, while the intensity values of almost all pixels of the Odd (Figure 5e) and Dbl components (Figure 5i) showed an increasing trend. This suggests that the absence of POAC will lead to an overestimation of Vol and an underestimation of Odd and Dbl. In addition, the correction of the ESAC and AVEC stages did not cause the pixel intensity value of the polarization decomposition component to undergo an overall increase or decrease, but there was an increase or decrease in the pixel value in each component. This is because the correction factors of ESAC and AVEC stages depend on both the value of the local incidence angle of each pixel and its relationship with the radar incidence angle, which also indicates that the correction effects of ESAC and AVEC stages are closely related to the ground conditions in the study area. However, in the ESAC phase, most of the pixels that cause a large change in the intensity value of the pixel (away from the 1:1 line) were decreasing, while in the AVEC phase, the situation was the opposite. This indicates that going without ESAC will lead to a more serious overestimation of some pixels, while going without AVEC will lead to a more serious underestimation of some pixels. Moreover, the effect of RTC on the three polarization decomposition components cannot be described by a simple law of increase and decrease because the correction effect of RTC is the superposition of the three stages (Figure 5d,h,l).

In addition, in order to comprehensively evaluate the universality and robustness of the effect of RTC on polarization decomposition, we used the Yamaguchi three-component (YAM3) method for analysis and map** (Figure A1, Appendix A). On the whole, the trend of the scatter density plots of the two polarization decomposition methods was the same in each stage except for slight differences caused by the different decomposition algorithms. This shows that the effect of RTC on the components of polarization decomposition is relatively robust in different polarization decomposition methods.

4.3. The Impact of Radiometric Terrain Correction on Regression Model Performance

In order to analyze the impact of RTC on the performance of machine learning models and select the optimal regression model, we conducted feature selection (Figure A2, Appendix A) and model training (

Table A2. The average value of the corresponding evaluation index of the optimal performance obtained by repeating 10 times: two-step feature selection was carried out using the radiometric terrain correction (RTC) and non-RTC (NRTC) datasets, and then 10 regression models and the Optuna algorithm were used to obtain the corresponding evaluation indicators of optimal performance.

Model	R2_RTC	R2_NRTC	rRMSE_RTC (%)	rRMSE_NRTC (%)
BysRidge	0.86334	0.69453	18.11756	22.89453
ARD	0.86213	0.71534	18.31867	22.57860
Lasso	0.85586	0.68453	18.13603	23.37506
ElasticNet	0.84783	0.68124	18.51316	22.56761
Ridge	0.84704	0.64753	18.74611	23.45273
CatBoost	0.82521	0.73013	18.14501	22.43760
XGBoost	0.82202	0.74273	19.01550	21.37638
AdaBoost	0.81934	0.69574	18.37307	22.45734
ET	0.79859	0.70750	18.56786	22.89376
RF	0.78988	0.68545	18.53375	22.45703

, Appendix B) for 10 regression models based on the NRTC data and RTC data (25 July 2020), respectively. Then, we conducted a descending order analysis and graphed the results (Figure 6) based on the R² value of the training results of the RTC data.

Compared with training using the RTC data (Figure 6a), R² values obtained by different models based on training using the NRTC data are generally lower (0.14 lower on average), while rRMSE values are generally higher (4.2% higher on average), indicating that not implementing the RTC process will cause greater estimation error in forest AGB inversion. In addition, from the R² values of different models trained on the NRTC data and RTC data, it can be found that the RTC process generally improves the performance of linear models more than non-parametric models, indicating that the existence of terrain factors will cause the real relationship between SAR information and forest AGB to become more complex and difficult to describe by a simple linear relationship. Meanwhile, RTC can effectively reduce the signal distortion caused by terrain factors and make SAR data more closely reflect the real situation of forest biomass, thus reducing the complexity of data.

Figure 6b is a scatter plot of the measured forest AGB of 132 plots and the predicted value obtained using the BysRidge regression model (the largest R² value in 10 repetitions of Optuna hyperparameter optimizations, Equation (11)), which had the best performance on the data of this study, indicating that the model has a strong ability to explain the data trend (R² = 0.82). Its prediction error was relatively small (rRMSE = 18.06%), but there were still some problems of overestimation and underestimation. In addition, the multivariate model has a higher R² value (0.82) compared to the R² value of the univariate model (0.73 for the data dated 25 July 2020 in our previous study [40]), indicating that the multivariate model has a stronger advantage in retrieving forest AGB. Moreover, we evaluated the generalization ability of the model (Equation (11)) with respect to the PolSAR data from different dates (11 July and 8 August 2020), and created a scatter plot of the measured values and the predicted values of forest AGB calculated based on the model (Figure A3). Finally, the forest AGB in the whole range of the study area was predicted and mapped (Figure 6c) based on the trained model.

$$\begin{array}{c}AGB(db)=56.4598+3.5369\times {\sigma }_{HV}^0\_db+1.5128\times BMI\_db+3.5913\times V3R\_2\_db\\ +23.2417\times H+2.0974\times mRFDI+2.7412\times F3V\_db-6.1481\times RVI\\ -2.9493\times Span\_db-0.5524\times R\_VV/HV\_db-26.397\times A3V\end{array}$$

(11)

5. Discussion

5.1. The Significance of Radiometric Terrain Correction

As an important part of the terrestrial ecosystem carbon pool [88], the accurate estimation of forest AGB is of great significance for quantifying the value of forest carbon storage and assessing the dynamics of the global carbon cycle, but the interference of topographic factors cannot be ignored [89]. This study found that the effective implementation of the RTC process can significantly suppress the interference of terrain factors with PolSAR data, thereby improving the quality and application potential of PolSAR data. This is specifically reflected in the following points. On the one hand, compared to the R-value between forest AGB and SAR features extracted from the NRTC data, the average R-value of 75 original features extracted from the RTC data increased by 0.26, and the maximum improvement was 0.48 (${\sigma }_{HV}^0$). On the other hand, RTC is highly reasonable and robust at correcting components obtained by different polarization decomposition methods (FRE3 and YAM3 are examples). Moreover, results consistent with those of other researchers [90,91,92] were obtained during the POAC stage: the Vol component decreased, while the Odd and Dbl components increased. However, the correction effect on each ground unit during the ESAC and AVEC stages was closely related to the ground conditions in the study area and did not exhibit a monotonic trend. Finally, the fitting performance of 10 regression models based on RTC data was generally superior to that based on NRTC data (R² increased by 0.14 and rRMSE decreased by 4.2% on average). These results indicate that RTC is effective and robust for improving the accuracy of PolSAR features and the reliability of forest AGB inversion. This is because RTC is a set of effective step-by-step terrain correction processes designed by in-depth analysis of the SAR physical mechanism, which accurately corrects the interference caused by terrain factors in terms of the scattering mechanism (POAC and AVEC) and intensity information (ESAC). The original PolSAR features after RTC processing carry more accurate polarization scattering information, thus improving the sensitivity of AGB inversion. In addition, if a suitable RTC process is created in the application of other SAR technologies to eliminate the interference of terrain factors, the application effect will be improved to a certain extent.

5.2. The Overcorrection of Radiometric Terrain Correction

This study also found that, without RTC processing, the R-value of the derived features extracted based on ratio arithmetic is significantly higher than that of the original features [86], but the original features showed significant improvement after RTC. This indicates that ratio-derived features can be prioritized for SAR data application when RTC is not feasible [93,94]. However, compared with the derived features without RTC, the average R-value of the 35 ratio-related PolSAR-derived features after RTC decreased by 0.05, indicating that there are some overcorrection issues in the RTC process that need to be improved. In addition, since the ratio operation can eliminate multiplicative noise, there is a greater likelihood of overcorrection due to POAC and AVEC. The extraction of POA needs to distinguish between the offset caused by forest targets and terrain factors, and only correct the offset caused by terrain, while the correction factor in the AVEC stage only considers the local incidence angle and radar incidence angle, without introducing factors such as surface roughness, terrain slope, and complex dielectric coefficient that are related to the backscatter coefficient in the AVEC correction process. Moreover, there is potential space for further improvement of the ESAC method. Therefore, exploring how to effectively solve the overcorrection problems caused by the above factors in the RTC process is a valuable future research direction.

5.3. Performance Comparison of Different Regression Models

Improving dataset quality (such as RTC) has a greater impact on improving AGB estimation performance than selecting regression algorithms [95]. However, after the significant improvement in data quality, the selection of machine learning regression models becomes particularly crucial. Selecting an appropriate regression model can ensure that the algorithm can fully utilize the features in the data, capture the relationships between variables, and avoid overfitting or underfitting, thereby further improving the prediction accuracy and generalization ability of the model. This study compared 10 different regression models and found that the linear model has superior complexity (feature dimension) and interpretability compared to the non-parametric model. This is because non-parametric models do not preset the specific form of the model but adjust flexibly based on the data. However, it is precisely because of their higher flexibility that these models become too complex. In the feature selection stage, the variance inflation factor (VIF) values of the features of different linear regression models are all less than 10 (Figure A2d), indicating that regularization methods (L1 or L2) can effectively alleviate the feature collinearity problem in multivariate linear models. However, the VIF value of many features is greater than 5, indicating that there is still a moderate degree of collinearity problem in the feature set. In addition, it is to be expected that there is some degree of collinearity in the features selected by the non-parametric model (Figure A2e), but the collinearity problem usually has less effect on the non-parametric model.

In addition, the R² and rRMSE indicators of the non-parametric models based on the RTC data in this study are usually worse than those of the linear models. This may be because the target variable of this study is the forest AGB value in decibels, which has a better linear relationship with PolSAR features (higher accuracy) [7,58,96]. However, non-parametric models based on the NRTC data in this study typically have better R² and rRMSE metrics than linear models. This is because terrain factors cause PolSAR features to have a more obviously nonlinear relationship with forest AGB. Moreover, the setting of hyperparameter spaces for different models directly affects the performance of the model, while the setting of hyperparameter spaces for non-parametric models is more complex, which may be one reason why non-parametric models performed slightly worse on the data in this study. Furthermore, the comparison between parametric and non-parametric models may be unfair because there is not a large amount of data to train the non-parametric model, resulting in a feature selection method that is not ideal. Therefore, k-Nearest Neighbor is a regression algorithm worth trying because it is more suitable for small datasets than random forest. In addition, deep learning models (such as Convolutional Neural Networks (CNN) and Deep Neural Networks (DNN)) have certain limitations (such as a large amount of data being required for training, the consumption of computing resources being large, and the model interpretation being poor) compared to traditional machine learning models. However, because of their strong feature learning ability, high accuracy, and excellent generalization ability, they have shown great potential for solving complex problems, and so are worth putting into use in scenarios with sufficient resources.

5.4. Potential Limitations

In this study, only PolSAR features were applied in AGB inversion. However, it is necessary to include more actual parameters that reflect the geographical environment and physical conditions when constructing a more comprehensive and accurate AGB model [44]. For example, climatic factors (temperature, rainfall, humidity, etc.) may directly affect the state of the land surface and vegetation, thereby altering the reflection characteristics of SAR signals [97]. Moreover, although RTC is carried out in the pre-processing of PolSAR data, it may also be worth taking slope, aspect, and local incidence angle into account in model training. In addition, more polarization parameters also contribute to the application of SAR data [67,98,99]. Through the comprehensive analysis of these multivariate parameters, we can expect a significant improvement in the predictive power of the model and a deeper understanding of the object of study.

In this study, data from only 132 plots were used to build and train machine learning models due to human constraints and real-world conditions. The relatively small forest AGB range (36.39–166.58 t/ha) of the real test sample may not be suitable for the forest in the study area (the maximum value is 254.6 t/ha). However, adding more plot data through unmanned aerial vehicle-light detecting and ranging (UAV-LiDAR) systems can significantly improve the generalization and reliability of the model, as larger datasets more accurately reflect the ecological diversity and complexity of the study area, hel** the model capture more patterns and trends [5]. Furthermore, the control of data quality is crucial. It is necessary to strictly control field measurement errors, AGB calculation errors, and matching errors between sample location and remote sensing image data to reduce the uncertainty caused by human factors, so as to comprehensively improve the overall quality of data and the accuracy of the model prediction [100].

In addition, this study did not conduct differentiated modeling for different forest types. However, the AGB modeling method based on forest type is a forest AGB modeling strategy that has attracted much attention in recent years. Because it is based on the basic characteristics of forest ecosystem diversity and heterogeneity, which is helpful to reveal the biomass distribution law and the fine modeling of forest AGB, so it has significant advantages in terms of the accuracy and practicability of the AGB inversion model [101,102].

6. Conclusions

Considering that the impact of RTC on the correctness of feature extraction and regression model performance needed to be further evaluated, this study took PolSAR data provided by ALOS-2 as the data source, and used 10 machine learning regression models and the Optuna hyperparameter optimization algorithm to verify the effectiveness and robustness of RTC in improving the correctness of the PolSAR feature set and the accuracy of forest AGB inversion. The detailed results are as follows: The RTC process effectively and robustly improved the correlation between the PolSAR features and the measured forest AGB (with an average increase of 0.26 in R-values) and the performance of the regression models (with an average increase of 0.14 in R² values and an average decrease of 4.2% in rRMSE). The RTC process also showed prime robustness and rationality for the correction of different polarization decomposition components. The ratio operation of PolSAR features offset some multiplicative noises (including terrain), which indicates that there was a certain degree of overcorrection in the RTC process. Therefore, the ratio class-derived features were given priority when there was no condition for RTC; after RTC, σ⁰_HV_db still had the highest R-value of all PolSAR features and forest AGB decibels. In addition, the PolSAR features had a closer linear relationship with the forest AGB of decibel states. In this case, the linear model remained a powerful and practical choice due to its efficiency and stability. For example, the optimal regression model in this study was BysRidge (R² = 0.82, rRMSE = 18.06%).

However, this study treated all forest cover areas as the same type without differentiated modeling of different forest types, and confirmed that there is an overcorrection problem in the RTC process. More PolSAR features and actual parameters, as well as more accurate sample data, should be applied in machine learning inversion of forest parameters. Although this study had these limitations, it also provides insight for future research. It is expected to significantly improve the inversion accuracy of forest AGB by overcoming the above limitations, which can provide more reliable data support for the scientific management and protection of forest resources.

In summary, this study validated the effectiveness and robustness of RTC in improving the correct extraction of PolSAR features and accurately inverting forest AGB, and analyzed the potential mechanism of RTC on polarization decomposition parameters. The method used in this study has good robustness and universality. It can therefore provide assistance and reliable reference values for the inversion of forest AGB and related SAR applications based on machine learning models using PolSAR data in complex terrain areas.