Soil Organic Carbon Estimation and Transfer Framework in Agricultural Areas Based on Spatiotemporal Constraint Strategy Combined with Active and Passive Remote Sensing

Abstract

Mapping soil organic carbon (SOC) plays a crucial role in agricultural productivity and water management. This study discusses the potential of active and passive remote sensing for SOC estimation modeling in agricultural areas, incorporating synthetic aperture radar (SAR) data (L-band quad-polarization and C-band dual-polarization), multi-spectrum (MS) data, and brightness temperature (TB) data. The performance of five advanced machine learning regression (MLR) models for SOC modeling was assessed, focusing on spatial interpolation accuracy and cross-spatial transfer accuracy, using two field observation datasets for modeling and validation. Results indicate that the SOC estimation accuracy when using MS data alone is comparable to that of using TB data alone, and both perform slightly better than SAR data. Radar cross-polarization ratio index, microwave polarization difference index, shortwave infrared reflectance, and soil parameters (elevation and soil moisture) demonstrate high correlation with the measured SOC. Incorporating temporal features, as opposed to single-phase features, allows each regression model to reach its upper limit of SOC estimation accuracy. The spatial interpolation accuracy of each MLR algorithm is satisfactory, with the Gaussian process regression (GPR) model demonstrating optimal modeling performance. When SAR, MS, or TB data are used individually in modeling, the estimation errors (RMSE) for SOC are 0.637 g/kg, 0.492 g/kg, and 0.229 g/kg for the SMAPVEX12 sampling campaign, and 0.706 g/kg, 0.454 g/kg, and 0.474 g/kg for the SMAPVEX16-MB sampling campaign, respectively. After incorporating soil moisture and topographic factors, the above RMSEs for SOC are further reduced by 57.8%, 35.6%, and 3.5% for the SMAPVEX12, and by 18.4%, 8.8%, and 3.4% for the SMAPVEX16-MB, respectively. However, cross-spatial transfer accuracy of the regression models remains limited (RMSE = 0.866–1.043 g/kg and 0.995–1.679 g/kg for different data sources). To address this, this study reduces uncertainties in SOC cross-spatial transfer by introducing terrain factors sensitive to SOC (RMSE = 0.457–0.516 g/kg and 0.799–1.198 g/kg for different data sources). The proposed SOC estimation and transfer framework, based on active and passive remote sensing data, provides guidance for high-resolution regional-scale SOC mapping and applications.

1. Introduction

Soil organic carbon (SOC), as one of the largest carbon reservoirs in terrestrial ecosystems, plays a critical role in soil and ecosystem functions [1,2]. By absorbing and storing atmospheric carbon dioxide, it helps mitigate climate change [3]. SOC is a major component of soil organic matter (SOM), improving soil pore structure, enhancing water retention, and increasing nutrient availability, thereby promoting plant growth [4,5]. SOC provides energy for soil microorganisms, stimulating microbial activity and biodiversity, which are essential for nutrient cycling [6]. Moreover, SOC helps stabilize soil structure, reduces erosion, and increases resistance to soil degradation [7,8]. Overall, SOC is not only vital for agricultural productivity but also plays a key role in global carbon cycling and climate regulation [3,9].

In situ field sampling and laboratory analysis of SOC often achieve high measurement accuracy, but these methods are limited to local sampling points and cannot adequately capture large-scale distribution trends [10]. Interpolation algorithms are generally used for mapping at regional scales, but they introduce considerable uncertainty [11,12]. In recent years, significant progress has been made in estimating SOC using remote sensing techniques [1,13,14,15]. However, challenges related to high spatial resolution and high-precision estimation still remain. For instance, the spatial resolution of currently available SOC gridded products is limited to 250 m or 1000 m [13], although resolutions of up to 90 m can be achieved in specific regions [14]. Moreover, gridded products often lack timeliness, which is particularly critical for agricultural production and management [15].

SOC remote sensing estimation typically integrates multi-source data such as multi-spectrum (MS), synthetic aperture radar (SAR), and LiDAR to obtain more comprehensive multidimensional information [16,17,18,19]. For example, optical data, especially hyperspectral data, can capture soil spectral characteristics, SAR data are more sensitive to soil structure and moisture, and LiDAR is more responsive to soil texture [20]. Traditional models, such as linear regression and principal component analysis (PCA), have gradually been replaced by more advanced machine learning regression (MLR) algorithms, leading to improved SOC estimation accuracy [21,22,23]. With the advancement of remote sensing platforms and the availability of high spatiotemporal resolution data, such as Sentinel and Landsat series, researchers are now better equipped to monitor the spatiotemporal dynamics of SOC at regional and even global scales [24,25,26].

Although hyperspectral remote sensing can provide detailed soil spectral characteristics, thereby improving the estimation accuracy and robustness of SOC, the limited availability of hyperspectral platforms currently restricts its application potential [17,27]. It is necessary to explore the potential of other multi-source remote sensing data, such as MS, active microwave, and passive microwave data (Brightness temperature, TB), for SOC modeling, as these data are more readily accessible and can comprehensively reflect SOC across different information dimensions. Moreover, the cross-spatial transferability of SOC estimation models is crucial for their application in SOC mapping and monitoring. Enhancing the cross-spatial generalization capability of SOC estimation models to ensure they are applicable to different soil types and regions has become an important research trend [28].

This study mainly discusses the potential of several MLR algorithms, combined with multi-source data, for estimating SOC in agricultural areas with complex surface heterogeneity. It develops an SOC estimation and transfer framework that combines active and passive remote sensing data with spatiotemporal constraints and evaluates its accuracy and uncertainty using two field-measured datasets. The study addresses the following key issues:

(1)

What is the significance of SAR (L- and C-band) features, MS features, TB features, and soil parameters (elevation, slope, soil moisture, and soil roughness) for estimating SOC?

(2)

Can incorporating multi-source data enable the regression algorithms to achieve optimal SOC estimation accuracy? Can temporal features improve the estimation accuracy of SOC? Which regression algorithm is most suitable for SOC modeling?

(3)

How to improve the regional-scale spatial transferability of the SOC regression model?

2. Study Areas and Data

2.1. Study Areas

The study area is located in Manitoba, Canada (Figure 1a,b), characterized by a typical temperate continental climate, with hot, sunny summers and long, cold winters. The region is flat and open, serving as one of the primary agricultural areas, where crops such as canola, wheat, oats, soybeans, and corn are predominantly grown (Figure 1f). In situ field data for this study were obtained from the SMAPVEX12 and SMAPVEX16-MB datasets (Soil Moisture Active Passive Validation Experiments 2012 and 2016-Manitoba) [29,30]. The area’s land surface is covered with diverse vegetation, features highly heterogeneous soil textures, and exhibits significant variations in surface soil moisture (SSM), making it an ideal environment for thoroughly testing SOC estimation models.

2.2. In Situ Sampling Data

2.2.1. SMAPVEX12 Dataset

In June and July 2012, the SMAPVEX12 sampling campaign included measurements of SSM, soil texture types (STT), soil roughness (SR), and various indicators of vegetation growth. Each sampling plot contained 2 surface SOM sampling points (0–5 cm deep), with a total of 55 plots sampled in early June [4]. A total of 110 SOM samples were collected during the SMAPVEX12 sampling campaign. STT data captured the proportions of sand, silt, and clay ($f_{sand}$, $f_{silt}$, and $f_{clay}$), while SR data encompassed the root mean square height (RMSH) and correlation length (CL). Since soil parameters (SOM, STT, and SR) remained relatively stable throughout the observation period, these measurements were taken only once [4,31]. Soil samples were collected manually by digging, with each sample having a diameter and depth of 4.7 cm and 4.6 cm. Soil samples were then dried in an oven at 105 °C for 24 h to measure SSM. Additionally, the particle size was measured using the pipette method. The dried soil samples were further ground and passed through a 0.5 mm sieve, then heated at 375 °C for 16 h to obtain SOM [4,32]. SR data were obtained by digital processing after shooting with a pinboard and digital camera. All sampled data were quality controlled [4]. The sampled crops included canola, cereals (wheat and oats), soybeans, corn, and pasture. Additional details are available on the SMAPVEX12 website (https://smapvex12.espaceweb.usherbrooke.ca, accessed on 14 July 2014).

2.2.2. SMAPVEX16-MB Dataset

In June and July 2016, the SMAPVEX16-MB sampling campaign was carried out across 50 farmland plots. The sampling method used was similar to that of the SMAPVEX12 dataset. Each plot contained 5–6 surface SOM sampling points, covering a total of 50 plots sampled in early June [29]. A total of 252 SOM samples were collected during the SMAPVEX16-MB sampling campaign. All sampled data were quality controlled. Vegetation growth conditions were also recorded, including crops such as canola, wheat, oats, soybeans, and corn. Further details can be found on the SMAPVEX16-MB website (https://smapvex16-mb.espaceweb.usherbrooke.ca, accessed on 10 July 2021).

Whether in SMAPVEX12 or SMAPVEX16-MB, there are significant differences in soil’s physical and geometric parameters between different plots (Figure 2). Since both sampling campaigns collected SOM, a coefficient conversion was performed to obtain SOC, as follows [33]: SOC = SOM/1.724.

2.3. Remote Sensing Data

2.3.1. Microwave Remote Sensing Data

During June to July 2012, a total of 13 phases L-band quad-polarization (quad-pol) SAR images were acquired using the uninhabited aerial vehicle synthetic aperture radar (UAVSAR) platform (Figure 1d), with ground-based synchronous observational sampling conducted for 12 phases. The L-band quad-pol SAR data consisted of horizontal–horizontal (HH), horizontal–vertical (HV), vertical–horizontal (VH), and vertical–vertical (VV) polarization modes (

Table 1. Remote sensing features used in the regression models.

Sampling Campaign	Feature Types	Feature Descriptions
SAMPVEX12	L-band SAR features	${\sigma }_{\mathrm{H}\mathrm{H}}$$, {\sigma }_{\mathrm{H}\mathrm{V} \mathrm{o}\mathrm{r} \mathrm{V}\mathrm{H}}$$, {\sigma }_{\mathrm{V}\mathrm{V}}$; $\mathrm{C}\mathrm{P}{\mathrm{R}}_{\mathrm{H}}^{\mathrm{L}\mathrm{-}\mathrm{b}\mathrm{a}\mathrm{n}\mathrm{d}}={\sigma }_{\mathrm{H}\mathrm{V}}/{\sigma }_{\mathrm{H}\mathrm{H}}$$, \mathrm{C}\mathrm{P}{\mathrm{R}}_{\mathrm{V}}^{\mathrm{L}\mathrm{-}\mathrm{b}\mathrm{a}\mathrm{n}\mathrm{d}}={\sigma }_{\mathrm{V}\mathrm{H}}/{\sigma }_{\mathrm{V}\mathrm{V}}$$, \mathrm{C}\mathrm{o}\mathrm{P}{\mathrm{R}}_{\mathrm{H}-\mathrm{V}}^{\mathrm{L}\mathrm{-}\mathrm{b}\mathrm{a}\mathrm{n}\mathrm{d}}={\sigma }_{\mathrm{H}\mathrm{H}}/{\sigma }_{\mathrm{V}\mathrm{V}}$, RIA [34,35]
	L-band TB features	$\mathrm{T}{\mathrm{B}}_{\mathrm{H}}$$, \mathrm{T}{\mathrm{B}}_{\mathrm{V}}$$, \mathrm{M}\mathrm{P}\mathrm{D}\mathrm{I}={\displaystyle \frac{\mathrm{T}{\mathrm{B}}_{\mathrm{V}}-\mathrm{T}{\mathrm{B}}_{\mathrm{H}}}{\mathrm{T}{\mathrm{B}}_{\mathrm{V}}+\mathrm{T}{\mathrm{B}}_{\mathrm{H}}}}$ [36]
	Optical features (SPOT4/5)	${\rho }_{\mathrm{G}\mathrm{r}\mathrm{e}\mathrm{e}\mathrm{n}}$$, {\rho }_{\mathrm{R}\mathrm{e}\mathrm{d}}$$, {\rho }_{\mathrm{N}\mathrm{I}\mathrm{R}}$$, {\rho }_{\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}}$; $\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}={\displaystyle \frac{{\rho }_{\mathrm{N}\mathrm{I}\mathrm{R}}-{\rho }_{\mathrm{R}\mathrm{e}\mathrm{d}}}{{\rho }_{\mathrm{N}\mathrm{I}\mathrm{R}}+{\rho }_{\mathrm{R}\mathrm{e}\mathrm{d}}}}$ [37]$, \mathrm{N}\mathrm{D}\mathrm{M}\mathrm{I}={\displaystyle \frac{{\rho }_{\mathrm{N}\mathrm{I}\mathrm{R}}-{\rho }_{\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}}}{{\rho }_{\mathrm{N}\mathrm{I}\mathrm{R}}+{\rho }_{\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}}}}$ [38]$, \mathrm{N}\mathrm{I}\mathrm{R}\mathrm{v}=\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}\times \mathrm{N}\mathrm{I}\mathrm{R}$ [39]
SAMPVEX16-MB	C-band SAR features	${\sigma }_{\mathrm{V}\mathrm{H}}$$, {\sigma }_{\mathrm{V}\mathrm{V}}$$, \mathrm{C}\mathrm{P}{\mathrm{R}}_{\mathrm{V}}^{\mathrm{C}\mathrm{-}\mathrm{b}\mathrm{a}\mathrm{n}\mathrm{d}}={\sigma }_{\mathrm{V}\mathrm{H}}/{\sigma }_{\mathrm{V}\mathrm{V}}$, RIA
	L-band TB features	$\mathrm{T}{\mathrm{B}}_{\mathrm{H}}$$, \mathrm{T}{\mathrm{B}}_{\mathrm{V}}$, MPDI
	Optical features (Sentinel-2A)	${\rho }_{\mathrm{B}\mathrm{l}\mathrm{u}\mathrm{e}}$$, {\rho }_{\mathrm{G}\mathrm{r}\mathrm{e}\mathrm{e}\mathrm{n}}$$, {\rho }_{\mathrm{R}\mathrm{e}\mathrm{d}}$$, {\rho }_{\mathrm{R}\mathrm{e}\mathrm{d} \mathrm{e}\mathrm{d}\mathrm{g}\mathrm{e}1}$$, {\rho }_{\mathrm{R}\mathrm{e}\mathrm{d} \mathrm{e}\mathrm{d}\mathrm{g}\mathrm{e}2}$$, {\rho }_{\mathrm{R}\mathrm{e}\mathrm{d} \mathrm{e}\mathrm{d}\mathrm{g}\mathrm{e}3}$$, {\rho }_{\mathrm{N}\mathrm{I}\mathrm{R}}$$, {\rho }_{\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}1}$$, {\rho }_{\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}2}$; $\mathrm{NDVI}, \mathrm{NDMI}, \mathrm{NIRv}, \mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{\mathrm{R}\mathrm{e}\mathrm{d} \mathrm{e}\mathrm{d}\mathrm{g}\mathrm{e}}={\displaystyle \frac{{\rho }_{\mathrm{N}\mathrm{I}\mathrm{R}}-{\rho }_{\mathrm{R}\mathrm{e}\mathrm{d} \mathrm{e}\mathrm{d}\mathrm{g}\mathrm{e}}}{{\rho }_{\mathrm{N}\mathrm{I}\mathrm{R}}+{\rho }_{\mathrm{R}\mathrm{e}\mathrm{d} \mathrm{e}\mathrm{d}\mathrm{g}\mathrm{e}}}}$$, \mathrm{N}\mathrm{D}\mathrm{S}\mathrm{I}={\displaystyle \frac{{\rho }_{\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}1}-{\rho }_{\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}2}}{{\rho }_{\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}1}+{\rho }_{\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}2}}}$ [40,41,42]

).

The products used were the multi-look ground project complex data, derived from the original single look complex (SLC) data with a spatial resolution of 1.66 m × 0.8 m, resulting in a final spatial resolution of 4 m × 6 m. The radar incidence angles (RIA) range from 20.0° to 65.0° (Figure 1c). The original backscattering coefficient data were normalized to 40° using histogram matching method [43]. The backscattering coefficient images were processed by mean filtering (5 × 5 window) and were resampled to 20–30 m spatial resolution using the nearest neighbor method to match the optical remote sensing data (Figure 3). Due to the low flight altitude of the UAVSAR platform, RIAs of different plots differ greatly, which will bring great uncertainty to the SOC modeling.

Sentinel-1A Interferometric Wide Swath (IW) single look complex (SLC) products during April to October 2016 (non-frozen stage) were chosen as the C-band dual-polarization (dual-pol) SAR data sources for analyzing backscattering coefficients in the study area. The preprocessing method is visually depicted in Figure 3. This study acquired dual-pol backscattering coefficients, which encompass both VH and VV polarization (Table 1). The RIAs were extracted from header files. To reduce the impact of RIA on the SOC models, the dual-pol backscattering coefficient images were processed using second-order cosine normalization and uniformly adjusted to 40° [44,45].

2.3.2. Optical Remote Sensing Data

The SPOT-4 and SPOT-5 Level-1A data were used as optical remote sensing data sources in 2012. After radiometric calibration, atmospheric correction, orthorectification, and geometric registration, four-band land surface reflectance (LSR) images with a spatial resolution of 10–20 m were obtained (Figure 1e), including two visible light bands (VIS; ${\rho }_{Green}$ and ${\rho }_{\mathit{Red}}$), one near-infrared band (${\rho }_{NIR}$), and one short-wave infrared band (${\rho }_{SWIR}$). The Sentinel-2A Level-2A data were selected as optical data sources in 2016. After preprocessing, the LSR images were obtained with 10–20 m spatial resolution, including three VIS bands (${\rho }_{Blue}$, ${\rho }_{Green}$, ${\rho }_{\mathit{Red}}$), three red edge bands (${\rho }_{red edge1}$, ${\rho }_{red edge2}$, ${\rho }_{red edge3}$), two NIR bands (${\rho }_{NIR}$ and ${\rho }_{nNIR}$), and two SWIR bands (${\rho }_{SWIR1}$ and ${\rho }_{SWIR2}$).

2.3.3. Brightness Temperature Data and Preprocessing

The TB data were derived from the passive and active L-band system (PLAS; 1.413 GHz), which was mounted on a fixed-wing aircraft. The imaging incident angle was set to 40°. The land surface vertical polarization TB (TB_V) and horizontal polarization TB (TB_H) in the study areas were mapped [29,30]. The spatial resolution of TB_V and TB_H was about 650 m (2012) and 450 m (2016), corresponding to the scale of a sampling plot. The microwave polarization difference index (MPDI) was derived from TB_V and TB_H (Table 1) [36]. Considering the joint modeling of multi-source remote sensing data, this study omitted some of the sampling plots that were not in the TB imaging ranges. The imaging information of the above data are shown in Tables S1 and S2 in the Supplementary Materials.

3. Methods

The technical process of this study is as follows (Figure 3): (1) The collected L- and C-band SAR data, TB remote sensing data, optical remote sensing data, and field observation data were preprocessed; (2) A correlation analysis was performed on temporal remote sensing features, the measured soil parameters (SSM and SR), terrain factors (DEM and slope), and the measured SOC; (3) SOC estimation models were constructed based on regression algorithms; (4) The SOC estimation and transfer accuracy of different MLR models were assessed under multiple strategy constraints.

3.1. Machine Learning Regression Algorithms

3.1.1. Ensemble Learning Regression Algorithms

Ensemble learning regression (ELR) leverages multiple regression models to enhance predictive accuracy and robustness. By aggregating predictions, ELR balances the strengths and weaknesses of individual models, reducing both bias and variance while improving overall performance. ELR includes techniques such as bagging (bootstrap aggregating) and boosting. Bagging generates subsets to train independent models, with their predictions averaged for the final output. Boosting, on the other hand, iteratively trains weak learners to correct errors, with the final prediction being a weighted average of the learners’ results. In this study, extreme random tree regression (ETR) and eXtreme Gradient Boosting regression (XGBR) were employed as the bagging and boosting algorithms, respectively, both of which have demonstrated great potential in the SOC modeling [46,47].

3.1.2. Support Vector Regression Algorithm

Support vector regression (SVR) is a regression method based on statistical learning theory [48,49]. Its core idea is to find an optimal regression hyperplane such that as many training data points as possible are close to this plane, while penalizing points that fall outside a certain error margin (commonly known as the ε-insensitive zone). This approach balances model complexity with prediction error. The SVR model focuses only on points far from the boundary of the optimal regression hyperplane (i.e., support vectors), while points within the ε zone do not affect the model. This mechanism simplifies the model and reduces the risk of overfitting. SVR model can handle nonlinear regression problems by using kernel functions (such as linear, radial basis function, and Gaussian kernels) to map data into a high-dimensional feature space, performing well in environments with small sample sizes and high-dimensional data [50].

3.1.3. Gaussian Process Regression Algorithm

Gaussian process regression (GPR) is a non-parametric algorithm used to model the underlying relationships in data. It relies on the principles of Gaussian processes, which enable it to model uncertainty in functions and make predictions for new data points based on previously observed data [51,52,53]. The GPR model utilizes flexible kernel functions that can adapt to various types of functional relationships between different data sources. This makes GPR a powerful regression method, capable of not only capturing complex relationships within the data but also providing insights into predictive uncertainties [53,54,55].

3.1.4. Neural Network Regression Algorithm

Deep neural network (DNN) is a feedforward neural network with multiple hidden layers [56]. DNN performs multiple layers of nonlinear transformations to automatically extract complex features from the input data. During training, DNN updates its weights gradually through the backpropagation algorithm to minimize prediction errors. Neurons in each layer map inputs to nonlinear outputs using activation functions, such as ReLU or Sigmoid, extracting higher-order features layer by layer. Deep neural network regression (DNNR) is also used in SOC modeling applications [57].

3.2. Modeling Features and Strategies

In the MLR models, the input features ($X_1, X_2, X_3,\dots$) include different remote sensing data (SAR, TB, or MS features), soil parameters, and terrain factors, and the output feature (Y) is the measured SOC. SAR features include different polarization backscattering coefficients and corresponding radar indices [34,35], TB features include dual-pol TB and MPDI, while MS features include reflectance in various bands and corresponding spectral indices [37,38,39]. For the SMAPVEX12 campaign, the remote sensing features used for modeling include L-band quad-pol SAR features, L-band TB features, and four-band MS features (Table 1). For the SMAPVEX16-MB campaign, the remote sensing features used for modeling include C-band dual-pol SAR features, L-band TB features, and ten-band MS features (Table 1). Additionally, the impact of introducing the measured SSM, the measured SR, and terrain factors (DEM and slope) derived from the Shuttle Radar Topography Mission (elevation and slope) on the SOC estimation model was also considered.

3.3. Model Construction and Optimization

The measured SOC data are randomly divided into two parts (Part A and Part B) according to different plots and sampling points. First, Part A (half of all samples) is used as the training set and Part B (the other half of all samples) as the test set. Then, Part B is used as the training set and Part A as the test set. The estimation results of the test sets from the two modeling processes were used for the final accuracy evaluation. The probability distribution curves of SOC between the training set and the test set are largely consistent (Figure 4a,b). This modeling partitioning strategy is primarily used to evaluate the estimation modeling performance of each regression algorithm for SOC. This estimation accuracy can be regarded as the spatial interpolation accuracy of the SOC regression model. Furthermore, for evaluating cross-spatial transfer accuracy, the measured SOC data are randomly divided into two parts according to different sampling plots. Half of the sampling plots are used as the training set, and the remaining sampling plots are used as the test set. Due to the differences between sampling plots, the probability distribution curves of SOC between the training set and the test set show significant variation.

The Bayesian optimization (BO) algorithm was used to perform hyperparameter (HP) optimization for each basic regression model, with the acquisition function evaluated based on expected improvement per second. This approach fully exploits the potential of data-driven models for SOC modeling [58]. It is important to note that HP optimization was conducted using a 10-fold cross-validation method within the training set, without involving any data from the test set. The optimization strategies for different regression models are described below (Table S3): For the ETR model, the optimizable HPs include maximum (max) features, max depth, minimum (min) samples leaf (split), min weight fraction leaf, and min impurity decrease [46]; for the XGBR model, the optimizable HPs include subsample, learning rate, Gamma value (complexity control), max depth, min child weight, and regularization parameters (Alpha and Lambda) [47]; for the SVR model, the optimizable HPs include Kernel function, Box constraint, Kernel scale, and Epsilon [50]; for the GPR model, the optimizable HPs include Sigma value, basis functions, kernel functions, and kernel scale [59]; for the DNNR model, the optimizable HPs include the number of neurons, maximum epochs, learning rates, and batch sizes, Lambda (L2 regularization), and activation functions [60].

3.4. Model Validation

The following indices were used to evaluate the estimation accuracy of the SOC models, including Pearson’s correlation coefficient (R), mean absolute error (MAE), and root mean square error (RMSE). The expressions are as follows:

$$R={\displaystyle \frac{Cov\left(X,Y\right)}{\sqrt{Var\left(X\right)Var\left(Y\right)}}}$$

(1)

$$MAE={\displaystyle \frac{1}{n}}\times {\displaystyle {\sum }_{i=1}^n\left|X_{obs,i}-X_{\mathit{model},i}\right|}$$

(2)

$$RMSE=\sqrt{{\displaystyle \frac{{\displaystyle {\sum }_{i=1}^n{\left(X_{obs,i}-X_{\mathit{model},i}\right)}^2}}{n}}}$$

(3)

where $\mathit{Cov}\left(X,Y\right)$ is the covariance of X and Y. $\mathit{Var}\left(X\right)$ and $\mathit{Var}\left(Y\right)$ are the variance of X and Y, respectively. n is the total number of verification data. $X_{obs,i}$ and $X_{model,i}$ are ith estimated SOC and the measured SOC, respectively.

4. Results and Discussions

4.1. Sensitivity Analysis of Multi-Source Remote Sensing Features and the Measured SOC

As shown in Figure 5a, for the SMAPVEX12 campaign, the correlation (R) between the measured SOC and L-band quad-pol SAR features in June is significantly better than in July. This is mainly because the crops in June were in the early growth stage with lower vegetation cover, allowing SAR data to penetrate better and capture the soil scattering signal [61]. However, the overall correlation is relatively low (R = −0.373–0.364), and most of the SAR features fail to pass the significance test across most observation phases (p > 0.05). Among the SAR features, the cross-polarization ratio ($\mathrm{C}\mathrm{P}{\mathrm{R}}_{\mathrm{H}}^{\mathrm{L}\mathrm{-}\mathrm{b}\mathrm{a}\mathrm{n}\mathrm{d}}$) is relatively more sensitive to SOC. The correlation between the measured SOC and L-band TB features is higher (R = −0.752–0.818), with MPDI being the most effective. In addition, there are good correlations between the measured SOC and MS data (R = −0.679–0.364). The correlation of single-band reflectance features is significantly better than that of derived spectral indices, and the SWIR bands show higher correlations. Moreover, there are strong correlations between the measured SOC and temporal SSM and DEM (R = 0.752, −0.744), while the correlations with slope and SR (RMSH and CL) are relatively low (p > 0.05).

For the SMAPVEX16-MB campaign, similar conclusions are drawn (Figure 5b). The correlation between the measured SOC and C-band dual-pol SAR features is relatively low (R = −0.353–0.490), with the cross-polarization ratio ($\mathrm{C}\mathrm{P}{\mathrm{R}}_{\mathrm{V}}^{\mathrm{C}\mathrm{-}\mathrm{b}\mathrm{a}\mathrm{n}\mathrm{d}}$) performing better than the single-polarization backscattering coefficient. Compared to L-band SAR data, the C-band SAR data have weaker penetration but better represented surface soil scattering signals. TB features still show high correlations (R = −0.512–0.656), with MPDI remaining the best indicator. The correlation between the measured SOC and MS features during the non-vegetated period (NDVI_mean±std = 0.093 ± 0.046) is significant (R = −0.800–−0.328; p < 0.001), with higher correlations in the two SWIR bands. Although the correlations are lower during the vegetation growth period, the two SWIR bands still show relatively high correlations (p < 0.001). Compared to SMAPVEX16-MB, the SMAPVEX12 sampling campaign experienced a complete “wet-to-dry” event, resulting in a greater temporal range of SSM variation. The measured SOC and temporal SSM at all observation phases exhibit high correlations (R = 0.660–0.875), as does the correlation with DEM (R = −0.772), while the correlations with slope and soil roughness (RMSH and CL) still remain low (p > 0.05). These results suggest that multi-source remote sensing features hold potential for SOC modeling.

4.2. Spatial Interpolation Accuracy Evaluation of SOC by Using Different Regression Algorithms

It is important to examine the upper and lower bounds and robustness of SOC estimation accuracy when using SAR data, TB data, MS data, soil parameters, and their combined modeling. This is crucial in practical applications, where simultaneous access to these data types may not always be possible. As shown in Figure 6a1, for the same regression algorithm (SMAPVEX12 sampling campaign), modeling with only L-band quad-pol SAR data (Group1) yields the lowest SOC estimation accuracy, followed by modeling with only MS data (Group2), whereas using only L-band TB data (Group3) achieves the highest accuracy. Combined modeling of SAR and MS data (Group4) significantly improves estimation accuracy over single data source models but remains lower than Group3. Among the algorithms, the GPR model performs best overall (Group4; R = 0.924, MAE = 0.336 g/kg, RMSE = 0.494 g/kg). For Group3, the metrics are R = 0.982, MAE = 0.178 g/kg, and RMSE = 0.229 g/kg. An evident correlation appears, that is, samples with higher SOC generally correspond to higher SSM, and the opposite also holds true. When the soil parameters (elevation and time-series SSM) are incorporated (Figure 6a2), SOC estimation accuracy markedly improves across all feature groups (R > 0.900), although performance differences between feature groups remain similar, with GPR still performing best. Figure 6a3 shows Group4’s results, without overestimations or underestimations (R = 0.959, MAE = 0.267 g/kg, RMSE = 0.365 g/kg). However, for some regression algorithms, Group4’s accuracy may be slightly lower than Group1 and Group2, indicating that a single remote sensing data source can suffice for SOC modeling when the soil parameters are included.

As illustrated in Figure 6b, compared to the SMAPVEX12 (SOC range: 2.0–5.5 g/kg), the SOC range in the SMAPVEX16-MB sampling campaign is broader (SOC range: 1.0–6.5 g/kg), making SOC modeling more challenging for the latter. Overall, using only C-band dual-pol SAR data for modeling (Group1) still results in the lowest SOC estimation accuracy for a given regression algorithm, while models using only MS data or TB data achieve similar accuracy. Compared to four-band reflectance SPOT4/5 data, ten-band reflectance Sentinel-2 data appears more suitable for SOC modeling. Similarly, when the elevation and time-series SSM factors are introduced, Group1’s accuracy improves significantly, with Group2 and Group3 showing smaller improvements (Figure 6b1,b2). Figure 6b3 shows that the GPR model performs best in modeling (R = 0.957, MAE = 0.335 g/kg, RMSE = 0.442 g/kg), and the scatter plots continue to reveal the close relationship between the measured SOC and the measured SSM.

Table S3 shows the hyperparameter values for each regression algorithm that achieved the best SOC estimation accuracy. It is worth noting that while excessive temporal features can lead to a “curse of dimensionality” and may reduce SOC estimation accuracy, this study finds that dimensionality reduction methods like PCA model do not enhance SOC estimation accuracy of regression models. On the contrary, introducing temporal features can yield an upper bound for SOC estimation accuracy. For the SMAPVEX12, modeling with temporal SAR features (12 phases; Table S1) reduces estimation errors by 5.1–20.4% compared to single phase SAR features, as it does for optical and TB features (4.2–10.8% and 3.3–8.6%). The same pattern holds for the SMAPVEX16-MB (26 phases; Table S2), where errors decrease by 7.5–22.9%, 5.5–13.3%, and 4.2–9.8% for each feature groups, respectively. Additionally, this study finds that including spectral indices derived from MS reflectance (Table 1) does not effectively improve SOC estimation accuracy and may even have a slight negative impact. Compared to using only reflectance features, when additional spectral indices are incorporated into the modeling, the estimation errors (RMSE) of the various regression algorithms increased by approximately 3.7–4.3%. Hence, original reflectance features are sufficient for meeting SOC regression modeling requirements.

4.3. Spatial Transfer Accuracy Evaluation of SOC by Using Different Regression Algorithms

Section 4.2 analyzes the interpolation performance of various regression algorithms for SOC modeling. However, it is also essential to investigate the cross-spatial transferability of SOC models, which is particularly crucial for SOC mapping. This analysis is conducted by using certain plots as the training set for SOC modeling, while the remaining plots serve as the test set. As shown in Figure 7a, for the SMAPVEX12 sampling campaign with the GPR model, the SOC accuracy in the validation set is acceptable when L-band quad-pol SAR data are used for modeling alone (R = 0.895, MAE = 0.376 g/kg, RMSE = 0.576 g/kg). However, the transfer accuracy is considerably lower (R = 0.607, MAE = 0.806 g/kg, RMSE = 1.043 g/kg). The same trend is observed when MS data are used, where validation accuracy (RMSE = 0.523 g/kg) is significantly better than transfer accuracy (RMSE = 0.866 g/kg). Although the inclusion of the terrain factor (DEM) markedly improves validation accuracy, there remains a considerable gap between validation and transfer accuracy (RMSE = 0.221 g/kg → 0.516 g/kg; RMSE = 0.199 g/kg → 0.457 g/kg; Figure 7a1–a4). Figure 7a5,a6 show the SOC validation accuracy (10-fold cross-validation accuracy) and test set accuracy (spatial transfer accuracy) of different regression algorithms with different remote sensing data for modeling. It is evident that, without the terrain factor, validation accuracy substantially exceeds transfer accuracy. The inclusion of the terrain factor significantly narrows this accuracy gap. Similar conclusions are reached for the SMAPVEX16-MB sampling campaign (Figure 7b). When C-band dual-pol SAR data are used alone, the SOC validation accuracy (RMSE = 0.616 g/kg) outperforms the lower transfer accuracy (RMSE = 1.420 g/kg). The trend is similar with MS data (RMSE = 0.423g/kg and 0.995 g/kg). Introducing the terrain factor substantially improves both validation and transfer accuracy (RMSE = 0.496 g/kg → 1.198 g/kg; RMSE = 0.438 g/kg → 0.799 g/kg; Figure 7b1–b4), although transfer results shows notable overestimation at low values and underestimation at high values, which is consistent across different regression algorithms (Figure 7b5,b6). Figure 7 also illustrates the potential of TB data in modeling the spatial transfer of SOC. However, its relatively low spatial resolution blurs the heterogeneity of SOC within the plots, leading to the illusion of low RMSEs between the estimated and observed values (with relatively low correlation).

This indicates that regression models trained on limited measured data struggle to meet the practical accuracy requirements for large-scale SOC mapping. When there is a substantial difference between the actual SOC distribution at the regional scale and the SOC at sampling points, the cross-spatial transfer accuracy of the regression models can decrease significantly, with transfer results potentially becoming unreliable. In other words, the SOC mapping performance of regression models depends on the reasonable distribution of sampling points [13,14]. The study area features relatively flat terrain with an elevation change in less than 30 m, characteristic of the Red River Valley [4]. Sandy soils in the west are slightly elevated by a moraine, separating them from the heavy clay soils in the east, which accounts for most of the elevation difference. Loam soils predominate in the south near Carman, with high SOC content (2–5%), strong water-holding capacity, and fertility. This also explains the close relationship observed between elevation and the measured SOC. Therefore, the improvement in SOC spatial transfer accuracy by including terrain factors, as observed in this study, may be regionally specific and limited. The inclusion of the SSM factor might also improve SOC spatial transfer accuracy. However, without measured and auxiliary information, obtaining high-resolution (spatial and temporal) SSM raster products remains challenging and uncertain [62,63].

4.4. Comparison with Other Studies

Currently, global-scale SOC modeling and mapping applications typically use low spatial resolution Moderate Resolution Imaging Spectroradiometer (MODIS) data (250–1000 m) [13,15,27]. For regional fine-scale studies, researchers often rely on single optical data sources (such as Landsat or Sentinel-2) for SOC estimation. For example, Zhang et al. explored the potential of various MLR models to estimate SOC in the JiangHan Plain using Landsat-8 data. The results showed that temporal NDVI data-based joint modeling more accurately estimated SOC than single-phase data (RMSE = 3.718 g/kg), but the overall correlation was relatively low (R = 0.625). However, due to minimal topographical variation in the study area, the research concluded that the inclusion of terrain factors had little effect on SOC modeling [64]. Castaldi et al. assessed the potential and uncertainty of SOC estimation in Germany’s agricultural areas using Sentinel-2 data. This study also highlighted the irreplaceable importance of the SWIR bands in SOC estimation modeling. It further concluded that a 20 m spatial resolution is sufficient to capture SOC heterogeneity between different plots in agricultural areas [65]. Zhou et al. and Wang et al. also demonstrated that combining SAR (C- and L-band) and optical data leads to more accurate SOC estimation in agricultural regions compared to using a single data source [19,66].

However, few studies have explored the potential and value of multi-temporal, multi-source data (including passive and active microwave remote sensing, optical remote sensing, and soil parameters) in SOC estimation modeling for agricultural areas, and even fewer studies have focused on the spatial transferability performance of SOC regression models. This study finds that the SOC modeling performance using TB features obtained from airborne platforms is good; however, its relatively low spatial resolution (450–650 m) limits the fine-scale monitoring of SOC in highly heterogeneous agricultural areas. The even lower spatial resolution (<10 km) of TB data from satellite platforms further constrains its potential for SOC modeling [67]. Therefore, time-series open-access SAR and MS satellite data with relatively high spatial resolution (10–30 m) are ideal data sources and low-cost means for SOC modeling in agricultural areas. In addition, this study addresses the issue of poor spatial transfer accuracy of SOC regression algorithms, which makes it difficult to meet the requirements for large-scale mapping. This issue was improved by introducing terrain factors that are highly sensitive to SOC. However, this approach may have regional limitations. Future research still needs to explore methods for effectively improving SOC spatial transferability and mapping capabilities without relying on specific factors.

5. Conclusions

This study developed an SOC estimation and transfer framework for agricultural areas based on MLR algorithms. It explores the SOC modeling performance and potential of active and passive remote sensing data, including L-band quad-pol SAR data, C-band dual-pol SAR data, L-band TB data, and MS data. Validation and assessment were conducted using two field-observed datasets, leading to the following conclusions:

(1)

Sensitivity analysis shows a strong relationship between the measured SOC and temporal SSM in both dry and rainy seasons, indicating its potential to reflect spatial variations in regional SSM. Additionally, cross-polarization ratio (SAR feature), MPDI index (TB feature), and shortwave infrared reflectance (MS feature) demonstrate high temporal correlation with the measured SOC. Thus, using active microwave, passive microwave, or optical data alone holds potential for SOC modeling.

(2)

The SOC estimation accuracy achieved with MS data alone is comparable to that obtained with TB data alone, both performing slightly better than SAR data. Introducing temporal features can bring optimal SOC estimation accuracy across all regression algorithms. The spatial interpolation of each regression algorithm is satisfactory, with the GPR algorithm achieving the best SOC modeling performance (RMSE = 0.365 g/kg and 0.442 g/kg for SMAPVEX12 and SMAPVEX16-MB sampling campaigns).

(3)

The cross-spatial transfer accuracy of MLR algorithms remains limited (RMSE = 0.866–1.043 g/kg and 0.995–1.679 g/kg for different data sources). To reduce uncertainties in cross-spatial SOC transfer, this study incorporates terrain factors sensitive to regional-scale SOC (RMSE = 0.457–0.516 g/kg and 0.799–1.198 g/kg for different data sources). The SOC estimation and transfer framework proposed in this study provides valuable guidance for high-resolution, regional-scale SOC mapping and applications, with substantial application potential for open-access Sentinel and NISAR satellites.