Improved Inversion and Digital Mapping of Soil Organic Carbon Content by Combining Crop-Lush Period Vegetation Indices with Ensemble Learning: A Case Study in Liaoning, Northeast China

1. Introduction

Soil organic carbon (SOC) is a critical component of the global carbon cycle and an essential indicator of soil health, directly influencing soil quality, crop productivity, and broader ecosystem functions [1,2,3] . As a fundamental element of soil fertility, SOC directly affects crop yields, soil water-holding capacity, and nutrient cycling efficiency. Moreover, the stable storage of SOC plays a vital role in mitigating the increase in atmospheric CO₂ concentrations [4], and its temporal and spatial dynamics are closely linked to land use, agricultural management, and climate change[5]. Therefore, accurately quantifying the spatial distribution of SOC is of both scientific and practical significance for achieving sustainable agricultural development, optimizing land management policies, and advancing carbon neutrality goals. In agricultural ecosystems, accurately mapping SOC content at a regional scale is essential for assessing soil fertility in farmland, optimizing agricultural management practices, estimating carbon sequestration potential, and informing climate change mitigation strategies [6].

Traditional SOC monitoring relies primarily on field sampling combined with laboratory analysis [7]. While highly accurate at the local scale, this approach is time-consuming, labor-intensive, and unsuitable for large-scale, dynamic assessments. The development of remote sensing technology has provided an efficient means for SOC estimation and digital mapping at regional scales, rendering large-area SOC monitoring increasingly feasible [8,9]。SOC contains various functional groups whose chemical bonds influence spectral reflectance, particularly in red and near-infrared (NIR) regions [10]. Strong correlations between the SOC content and reflectance within the 600–800 nm spectral range have been widely documented [10,11,12], forming a solid theoretical foundation for SOC inversion using optical remote sensing data. Numerous case studies have further demonstrated the effectiveness and applicability of this approach in diverse agroecological environments [13] .

Accurate SOC inversion depends primarily on the following two aspects: (1) the construction of inversion models on the basis of the selection of predictive indicators and (2) the establishment of suitable inversion algorithm models[14,15].The existing inversion models generally incorporate two types of indicators—direct observation indicators and indirect influencing indicators. Direct observation indicators are mostly quantitative remote sensing indices derived from bare-soil imagery that use data-driven modelling methods, such as spectral mathematical transformations [16] and vegetation indices (VIs), such as the normalized difference vegetation index (NDVI), enhanced vegetation index (EVI), and clay index (CI) [14,15,16,17,18,19]. Indirect influencing indicators are typically derived from geographic information system (GIS)-based spatial data, including topographic and geomorphological parameters [14], precipitation [20],temperature[21] and other environmental variables [14,18,22,23]. A critical issue arises with direct indicators: many studies inappropriately incorporate VIs (e.g., the NDVI and EVI) during bare-soil periods when vegetation is scarce or absent [24,25,26,27,28,29]. VIs are inherently designed to reflect vegetation growth; thus, their inclusion in bare-soil models lacks ecological relevance—it introduces noise rather than meaningful information, despite the use of both SOC and VIs in the red and NIR bands.

This disconnect highlights a critical issue: vegetation itself is deeply intertwined with SOC dynamics. Vegetation drives organic matter input via litter and root exudates [24,30], whereas SOC fertility directly influences vegetation growth [31]. Recent studies have confirmed strong correlations between the SOC content and crop growth, with vegetation regulating microbial activity and SOC mineralization [32,33]. Thus, VIs from the crop-lush period (VIs_lush) could serve as ecologically meaningful indirect indicators of SOC, capturing the cumulative effect of SOC on vegetation productivity. In contrast to bare-soil period VIs (VIs_bare), VIs_lush reflect real-time vegetation status, which is inherently linked to SOC via soil–vegetation feedback. This ecological relevance suggests that VIs_lush may outperform VIs_bare in terms of SOC prediction. However, this hypothesis remains untested, and the literature lacks systematic comparisons of such indicator choices. Additionally, in terms of modelling approaches, existing algorithms can be broadly classified into mathematical models [18,34,35,36,37,38] and machine learning algorithms (MLAs). MLAs generally outperform traditional statistical methods in terms of accuracy, particularly when complex datasets with heterogeneous statistical distributions are used[39,40]. Among MLAs, ensemble learning models (ELMs), such as the random forest (RF), extreme gradient boosting (XGBoost), and adaptive boosting (AdaBoost) models, have yielded promising results for SOC content estimation [41,42,43,44]. However, their performance under different indicator strategies (e.g., VIs_lush vs. VIs_bare) is poorly understood.

According to the above discussion, vegetation cover is associated with increased organic matter inputs to the soil [24,30], whereas changes in soil properties can alter vegetation growth patterns. Therefore, VIs_lush are likely important indirect indicators of the SOC content [45]. We proposed that VIs_lush may serve as more effective predictors of SOC content than VIs_bare because of their closer ecological connections with vegetation productivity and organic carbon input processes. To test this hypothesis, Landsat 8 imagery from different phenological stages, representing both the bare-soil period (April) and the crop-lush period (August), was adopted. In addition, we integrated topographic, precipitation, and hydrological variables to construct a comprehensive set of predictive indicators. These factors were categorized into direct observation factors and indirect influencing factors, thereby improving the interpretability and predictive power of the inversion models. Liaoning Province, a major agricultural region in northeastern China, was selected as the study area to determine whether replacing VIs_bare with VIs_lush could increase the accuracy and stability of SOC inversion.

The specific objectives of this study are as follows: (1) to establish SOC inversion models by integrating both direct observation factors and indirect influencing factors; (2) to design and compare three modelling strategies (MSs), namely, combining other features with VIs_lush (MS-1); combining other features with VIs_bare (MS-2); and using other features without VIs (MS-3); (3) to evaluate the predictive performance of three ELMs—XGBoost, RF, and AdaBoost—under each MS and to apply the SHapley Additive exPlanations (SHAP) method to interpret model outputs and quantify the relative contributions of VIs_lush and VIs_bare; and (4) to identify the optimal model and MS and subsequently generate a spatially explicit SOC content distribution map of Liaoning Province, thereby revealing regional SOC patterns and providing a systematic reference for sustainable land and soil management.

2. Materials and Methods

2.1. Study Area and Soil Samples

Liaoning Province, which is located in the southern part of Northeast China, is among the key grain-producing provinces [46]. Geographically, it spans the area range from 38°43′N to 43°26′N and 118°53′E to 125°46′E, covering an area of approximately 148,000 km² and supporting a permanent population of approximately 43 million. The province exhibits a diverse topography, with elevations ranging from sea level to approximately 1,323 m. In general, the terrain is higher in the north and lower in the south, characterized by mountainous and hilly areas surrounding the central plains that gently slope towards the coast. The Liaoxi region (west of the Liao River) and the Liaodong region (east of the Liao River) are dominated by hills and mountains, with average elevations of approximately 800 and 500, respectively. The central area, referred to as the Liaohe Plain, has an average elevation of approximately 200 meters (Figure 1). This varied terrain significantly influences the climate, hydrology, and agricultural productivity of the province, suggesting its role as a vital agricultural base in China.

In Liaoning Province, mountainous regions occur in both the eastern and western areas. The eastern mountains comprise branches of the Chosan Range, such as the Badailing and Longgang Mountains, extending from North Korea and running north to south, with elevations ranging from 500 to 800 m. In the west, the Inner Mongolian Plateau gradually descends into the Liaohe Plain, forming a mountainous zone with elevations between 300 and 1,000 m. The land cover of Liaoning Province is composed of approximately 88,000 km² of mountainous terrain (59.5% of the total area), 48,000 km² of plains (32.4%), and 12,000 km² of water bodies and other land uses (8.1%). Cropland covers approximately 36,900 km², accounting for 38% of the total land area (Figure 1(b)). Liaoning Province has a mild climate, ranging from a temperate humid continental type in the south to a continental monsoon type in the west. The annual growing season precipitation varies from 450 to 1,200 mm, and the mean growing season temperatures range from 4.6°C to 10.3°C. The major soil types include brown soils, meadow soils, cinnamon soils, and paddy soils. The province’s primary crops are annuals such as corn and rice, with a typical bare-soil period occurring from early to late April before planting [3,44]. In this study, a total of 468 soil samples were collected from cropland fields across Liaoning Province to determine the SOC content (Figure 1(b)). Sampling was conducted in 2020, with each sample obtained at depths from 0–20 cm. Following collection, the samples were dried naturally and finely ground using a ball mill. The SOC content was then measured using an elemental analyzer (Thermo Fisher Scientific, USA) following the procedures outlined in [47] and [48]. The spatial distribution of the sampling sites is presented in Figure 1.

2.2. Data Acquisition and Treatment

The dataset used for model training included 468 SOC observations, which is relatively limited for machine learning applications. To enhance spatial representativeness and model robustness, ordinary kriging interpolation was applied within a 50-m radius, expanding the dataset to 2,799 sample points, as detailed in Table 1.

Remote sensing and GIS data were collected from both ​​direct observation​​ and ​​indirect influence perspectives​​. Remote sensing imagery from two key phenological periods in 2020 was employed: April (bare-soil period) for the direct observation inputs and August (crop-lush period) for the indirect influencing factors. The remote sensing data were primarily sourced from the Geo-Spatial Data Cloud of the National Geomatics Center of China (https://www.gscloud.cn) and the United States Geological Survey (USGS, https://earthexplorer.usgs.gov/), with a cloud cover threshold set to ≤20%. Specifically, Level-2 surface reflectance products from the Landsat 8 Operational Land Imager (OLI) were acquired for the selected periods. After acquisition, the imagery underwent preprocessing steps, including mosaicking and stitching. In cases where the imagery from 2020 exhibited excessive cloud contamination, supplemental images from the corresponding phenological period in 2021 were incorporated to ensure data completeness.

The processed images were then clipped to the administrative boundaries of Liaoning Province using vector data within the ENVI 5.6 software. The final composite images for the bare soil and peak vegetation periods are shown in Figure 2.

Elevation, slope, and surface runoff data were derived primarily from calculations based on digital elevation model (DEM) data. The DEM for Liaoning Province, with a spatial resolution of 30 m, was obtained from the National Geospatial Data Cloud (https://www.gscloud.cn). Land use data for the year 2020, with a spatial resolution of 30 metres, were sourced from Wuhan University (http://irsip.whu.edu.cn/resources/CLCD.php). Precipitation data, featuring a spatial resolution of 1 km, were acquired from the Resource and Environmental Science Data Center of the Chinese Academy of Sciences (https://www.resdc.cn/).

2.3. Construction of SOC Content Quantitative Inversion Model

Under the guidance of soil formation theory, a quantitative inversion model for the SOC content was developed on the basis of two categories of explanatory variables, namely, direct observation factors and indirect influencing factors (Figure 3). Direct observation factors were derived from spectral transformations and indices extracted during the bare-soil period, aiming to directly capture the spectral responses associated with SOC content levels. Indirect influencing factors were obtained from remote sensing imagery acquired during the crop-lush period, as well as from DEM data and climate data, with the objective of representing the environmental conditions that influence SOC formation, accumulation, and depletion.

To evaluate the model performance, three widely adopted ELMs—the RF, XGBoost, and AdaBoost models—were compared. The algorithm that achieved the best accuracy was selected to establish the final SOC content inversion model for croplands across Liaoning Province.

2.3.1. Direct Observation Factors

In accordance with the agricultural calendar of Liaoning Province, satellite imagery acquired during the bare-soil period (typically in April) was employed to extract direct observational factors. These factors included band reflectance [11,49]. Specifically, reflectance values from the red, NIR, shortwave infrared 1 (SWIR1), and shortwave infrared 2 (SWIR2) bands were extracted. In addition, three spectral indices were calculated: the normalized difference water index (NDWI) [30], the soil clay index (CI) [50,51], and the soil brightness index (BI) [17], which together characterize the spectral features associated with surface soil properties. These indices were computed formulas as follows:

$$NDWI={\displaystyle \frac{R_{green}-R_{nir}}{R_{green}+R_{nir}}}$$

(1)

$$CI={\displaystyle \frac{R_{SWIR1}}{R_{SWIR2}}}$$

(2)

$$BI=\sqrt{{R_{red}}^2+{R_{nir}}^2}$$

(3)

where, $R_{green}$, $R_{red}$, $R_{nir}$, $R_{SWIR1}$ and $R_{SWIR2}$ represent the reflectance values of the green, red, NIR, SWIR1 and SWIR2 bands, respectively.

The spatial distributions of the NDWI, CI, and BI values in Liaoning Province are shown in Figure A1(a), Figure A1(b) and Figure A1(c), respectively.

2.3.2. Indirect influencing factors

On the basis of the theoretical framework of soil formation models, we developed a set of indirect environmental factors by integrating VIs derived from remote sensing imagery during the crop-lush period (typically in August) [52,53], along with precipitation data [20,54], surface runoff features[20,55,56] and terrain attributes extracted from the DEM [57,58,59,60,61]. Three VIs were derived from Landsat 8 imagery acquired during August, namely, the NDVI, EVI, and the green chlorophyll index (GCI), and the calculation equations are provided below.

$$NDVI={\displaystyle \frac{R_{nir}-R_{red}}{R_{nir}+R_{red}}}$$

(4)

$$EVI=2.5\times {\displaystyle \frac{R_{nir}-R_{red}}{{1+R}_{nir}+6\times R_{red}-7.5\times R_{blue}}}$$

(5)

$$GCI={\displaystyle \frac{R_{nir}}{R_{green}}}-1$$

(6)

where, $R_{nir}$, $R_{red},$ $R_{blue}$​ and $R_{green}$ represent reflectance values of the respective bands.

The spatial distributions of VIs_lush (NDVI, EVI and GCI) across Liaoning Province are shown in Figure A2(a), Figure A2(b), and Figure A2(c), respectively. To assess the effectiveness of VIs during the crop-lush period for SOC prediction, a comparative analysis was performed using VIs from both the bare-soil period and the crop-lush period. The results for VIs_lush, which are shown in Figure A1(d), Figure A1(e), and Figure A1(f), provide a foundation for model selection and validate the contribution of VIs to SOC content inversion.

Surface runoff pathways were derived via the D8 flow direction algorithm applied to the DEM data, as shown in Figure A2(d). On the basis of these pathways, two additional hydrological features were extracted, namely, surface runoff density (SRD) and surface runoff buffer (SRB), as presented in Figure A2(e) and Figure (f), respectively.

Furthermore, the topographic wetness index (TWI) was calculated to quantify the potential for water accumulation on the basis of the terrain features. The TWI is widely employed as a terrain-derived proxy for surface soil moisture distribution and runoff potential, both of which are closely linked to SOC variability. The TWI was computed as follows:

$$TWI=\ln {\displaystyle \frac{CA}{\tan \theta }}$$

(7)

Where CA is the upslope contributing area (flow accumulation), and $\theta$ denotes the local slope angle. The spatial distribution of the TWI values in Liaoning Province is shown in Figure A1(h).

2.3.3. SOC Content Quantitative Inversion Model

Based on the analytical framework, a quantitative inversion model for SOC content in croplands across Liaoning Province was constructed. The model integrates both direct observation factors and indirect influencing factors, as summarized in Table 2. Corresponding GIS-based spatial distribution layers for each controlling variable are illustrated in Figure A1 and Figure A2.

2.4. Modelling Strategy

2.4.1. Ensemble Learning Algorithms

Ensemble learning increases predictive and classification performance levels through the integration of multiple models. By combining weak learners into a strong composite model, compared with single-model approaches, ELMs improve generalization and reduce the risk of overfitting [62,63]. The two main strategies in ELMs are bagging and boosting (Figure 4). Bagging employs parallel computation and bootstrap sampling to create diverse training subsets. Each subset is used to train an independent model, and the final prediction is obtained through averaging or voting, which reduces variance and enhances accuracy [64]. In contrast, boosting follows a sequential training process. Each model focuses on correcting the errors of its predecessor by adjusting sample weights, thereby reducing bias over iterations [65,66].

The RF model is a classic bagging algorithm that is renowned for its robustness and high performance across diverse datasets [67]. Boosting methods include the AdaBoost model, which adaptively adjusts sample weights to increase the prediction accuracy[66], and the XGBoost model, which employs gradient-based optimization techniques to enhance the model performance [68,69].

a) Random forest

The random forest model, developed by Leo Breiman in 1996, is an ensemble learning algorithm that employs the bagging strategy [64]. It was further enhanced through the introduction of the random subspace method by Tinkham in 1998 [67,70], which increases model diversity by randomly selecting subsets of features for each tree. Multiple decision trees are constructed [71], and final predictions are generated by averaging individual tree outputs:

$$\widehat{f}\left(x\right)={\displaystyle \frac{1}{M}}{\sum }_{t=1}^M{\widehat{f}}_t\left(x\right)$$

(8)

Where, M is the total number of decision trees, ${\widehat{f}}_t\left(x\right)$ is the prediction from the t-tℎ tree, and $\widehat{f}\left(x\right)$ is the final ensemble prediction.

The RF model is highly effective in feature selection, resilient to noise, and capable of parallel processing. However, it may underperform with biased datasets or when model parameters are not well tuned.

b) AdaBoost

The AdaBoost model provides an increased regression accuracy by iteratively training a series of weak learners. It adjusts the weights of training samples at each iteration, increasing the emphasis on the previous prediction error of its predecessor [72]. The final strong learner is a weighted combination of weak learners:

$$F(x)={\sum }_{t=1}^T{\alpha }_t{f}_t\left(x\right)$$

(9)

Where, $F\left(x\right)$denotes the final prediction, $f_t\left(x\right)$is the output from the $t$-th regression model, and ${\alpha }_t$represents the weight assigned based on its predictive accuracy. The weight of each weak learner is computed based on its average loss $L_t$ , using:

$${\alpha }_t={\displaystyle \frac{1}{2}}\ln ({\displaystyle \frac{1-L_t}{L_t}})$$

(10)

Where, $L_t$ is the regression error of the $t$-th weak learner. The sample weights are updated using:

$${w_i}^{(t+1)}={w_i}^{\left(t\right)}·\mathrm{e}\mathrm{x}\mathrm{p}({\alpha }_t·\left|y_i-f_t\left(x_i\right)\right|)$$

(11)

Where,${w_i}^{\left(t\right)}$is the weight of sample $i$ at iteration $t$, $y_i$ is the observed SOC content, and $f_t\left(x_i\right)$ is the model prediction result for sample $i$.

The AdaBoost model can capture complex nonlinear patterns, but its sensitivity to noisy observations and outliers remains a limitation in regression contexts.

c) XGBoost

The XGBoost model, developed by Chen and Guestrin [68], enhances traditional boosting techniques through second-order Taylor expansion of the loss function and the incorporation of regularization terms, which improves accuracy and reduces overfitting [63]:

$$L^{\left(t\right)}=\sum _{i=1}^n\left[l\left(y_i,{{\widehat{y}}_i}^{t-1}+f_t\left(x_i\right)\right)\right]+\Omega \left(f_t\right)$$

(12)

$$\Omega \left(f_t\right)=\alpha T+{\displaystyle \frac{1}{2}}\lambda {\sum }_{j=1}^T{\omega }_j^2$$

(13)

Where, $y_i$ is the true SOC content, ${{\widehat{y}}_i}^{\mathrm{t}-1}$ is the predicted value at the previous iteration, $f_t\left(x_i\right)$ is the prediction from the new regression tree added at $t$-th iteration , $l$ is a convex loss function, $\Omega \left(f_t\right)$ Is the regularization term that penalized model complexity, $T$ is the number of leaf nodes in the regression tree, $w_j$ is the score assigned to the $j$-th leaf, $\alpha$And $\lambda$ are regularization parameters that control tree complexity and weight magnitude, respectively.

In regression scenarios, the XGBoost model excels at modelling nonlinear relationships, mitigating overfitting, and managing high-dimensional input features efficiently.

2.4.2. SHapley Additive Explanations

To quantify the contribution of each input variable to the predicted SOC content, the SHAP framework was applied to the ELMs. The SHAP method originates from cooperative game theory, wherein it attributes the difference between a specific prediction and the average prediction to individual features [73]. On the basis of the foundational work of SHAP [74,75,76], this approach fairly allocates the model output among all features on the basis of their marginal contributions across all possible feature combinations.

Mathematically, the shapley value for feature $m$ is defined as:

$${\varnothing }_m\left(v\right)={\sum }_{S\subseteq N\left\{m\right\}}{\displaystyle \frac{\left|S\right|!\left(\left|N\right|-\left|S\right|-1\right)!}{\left|N\right|!}}(v(S\cup \{m\})-v(S))$$

(14)

Where, ${\mathrm{\varnothing }}_m\left(v\right)$ denotes the contribution of feature $m$, $N$ is the set of all features, $S$ is any subset of $N$ not containing $m$, $v(S\cup \{m\}$) and $v\left(S\right)$ represent the model output for the subsets with and without feature $m$, respectively. The term ${\displaystyle \frac{\left|S\right|!\left(\left|N\right|-\left|S\right|-1\right)!}{\left|N\right|!}}$ denotes the probability corresponding to various feature combinations. Therefore, the result obtained from the above expression represents the marginal contribution of feature $m$ to the final outputs.

When it comes to the total contribution of all features for each observation, The total prediction $g\left(x^{\mathrm{\text{'}}}\right)$ can be expressed as:

$$g\left(x^{\text{'}}\right)={\varnothing }_0+{\sum }_{m=1}^M{\varnothing }_m{z_m}^{\text{'}}$$

(15)

Where, ${\mathrm{\varnothing }}_0$ represents the model outputs without any features, $M$ is the number of all features, and ${z_m}^{\mathrm{\text{'}}}$ is the value of feature $m$ for the given observation.

2.4.3. Model Evaluation

To evaluate the ability of the three ELMs to predict the SOC content, we divided the SOC samples into a training dataset and a validation dataset at a 2:1 ratio. The final meta-model was trained using 10-fold cross-validation to tune and determine the optimal parameters of the model. The trained regression models were subsequently validated using a validation dataset, and their performance was evaluated through the Pearson correlation coefficient (R), coefficient of determination (R²), root mean square error (RMSE), ratio of performance to deviation (RPD) and ratio of performance to the interquartile distance (RPIQ) as follows:

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

(16)

$$R^2=1-{\displaystyle \frac{\sum _i{(\widehat{y_i}-y_i)}^2}{\sum _i{(\overline{y_i}-y_i)}^2}}$$

(17)

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

(18)

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

(19)

$$RPIQ={\displaystyle \frac{IQR}{RMSE}}$$

(20)

Where $n$ is the number of inversion units, $\widehat{y_i}$ is the predicted value of the i-th inversion unit, $y_i$ is the true value of that inversion unit, $\overline{y_i}$ is the mean of observed values, and $\overline{\widehat{y_i}}$ is the mean of predicted values, IQR is the inter-quartile range of the observed values, calculated as the difference between the 75th percentile (Q3​) and the 25th percentile (Q1), IQR = Q3-Q1.

Both R and R² reflect the goodness of fit of models; R is particularly useful for evaluating the linear association, whereas R² emphasizes the proportion of explained variance. The RMSE quantifies the average magnitude of prediction errors and is sensitive to large deviations due to the squared term. RPD is the ratio of the standard deviation of observed values to the RMSE, whereas the RPIQ relies on the interquartile range (IQR) instead of the standard deviation, making it more robust to outliers and skewed distributions. Smaller RMSE values, along with larger RPD and RPIQ values, indicate greater predictive ability of the model [77]. These metrics provide a comprehensive and complementary evaluation of model accuracy, robustness, and generalizability.

2.4.3. Three Modelling Strategies

In this study, three MSs were proposed on the basis of direct observation and indirect influencing factors (Table 3). These strategies were designed to evaluate and compare the predictive performance of VIs derived from different phenological stages as follows:

Modeling Strategy 1 (MS-1): VIs_lush are used along with other factors.

Modeling Strategy 2 (MS-2): VIs_bare are used in conjunction with the same set of additional factors.

Modeling Strategy 3 (MS-3): Modelling is conducted without VIs.

3. Results

3.1. Accuracy of the Different Modelling Strategies and Ensemble Learning Models

The predictive performance of the three ELMs under MS-1, MS-2, and MS-3 is shown in Figure 5 and Figure 6. A comparative analysis of these strategies indicates that the two VI-enhanced strategies (MS-1 and MS-2) provided a significantly enhanced model performance. Regardless of the ensemble algorithm used, MS-1 consistently provided the most accurate predictions. This finding demonstrates that incorporating vegetation information, especially from crop-lush periods, provides essential ecological insights and enhances the model's ability to capture spatial variations in SOC. Overall, the R² values of all three ELMs greatly increased with the addition of VIs.

Under MS-1, the XGBoost model achieved an R² value of 0.84, an RMSE of 2.22 g/kg, an RPD of 2.49 and an RPIQ of 3.25. In contrast, under MS-2, the performance decreased to an R² value of 0.79, an RMSE of 2.53 g/kg, an RPD of 2.18 and an RPIQ of 2.85. The accuracy further decreased under MS-3, with an R² value of 0.77, an RMSE of 2.63 g/kg, an RPD of 2.10 and an RPIQ of 2.74. These results show that compared with MS-3, the inclusion of VIs_bare under MS-2 increased the R² value of the XGBoost model by 2.60%, reduced the RMSE by 3.80%, increased the RPD by 3.81%, and enhanced RPIQ by 4.01%. The use of VIs_lush in MS-1 led to a 9.09% increase in R², a 15.59% decrease in RMSE, an 18.57% increase in the RPD, and an 18.61% increase in RPIQ. When MS-1 was compared with MS-2, the R² improved by 6.33%, the RMSE decreased by 4.00%, the RPD increased by 18.61%, and the RPIQ increased by 14.04%. These improvements highlight the high explanatory power of VIs, especially VIs_lush, in modelling SOC variability. With respect to the RF model, MS-1 yielded an R² value of 0.70, an RMSE of 3.04 g/kg, an RPD of 1.82, and an RPIQ of 2.37, which were noticeably better than those under MS-2 (R² = 0.67, RMSE = 3.17 g/kg, RPD = 1.74, and RPIQ = 2.27) and MS-3 (R² = 0.67, RMSE = 3.16 g/kg, RPD = 1.75, and RPIQ = 2.28). However, it is noteworthy that in the RF model, the results of MS-1 and MS-2 were quite close, with even slight decreases observed in some metrics (e.g., RMSE, RPD, and RPIQ). This suggests that using VIs_bare may not enhance SOC prediction in this case and might even introduce data redundancy. From MS-3 to MS-2, the R² remained unchanged, whereas from MS-2 to MS-1, the value increased by 4.48%. This further confirms that VIs_lush provide important information that is not captured by the terrain or soil texture variables alone. Although the AdaBoost model generally showed lower predictive ability, a similar pattern was observed. Its prediction results resembled those of the RF model, with little to no change in R², RMSE, RPD, or RPIQ from MS-3 to MS-2 but with slight improvements upon adding VIs_lush under MS-1. Specifically, R² increased from 0.33 (MS-3) to 0.34 (MS-1), representing a 3.03% improvement, the RMSE decreased from 4.53 g/kg to 4.51 g/kg (−0.44%), the RPD increased from 1.22 to 1.23 (+0.82%), and the RPIQ increased from 1.59 to 1.60 (+0.63%). These results again underscore the crucial role of VIs_lush in enhancing prediction accuracy, albeit with more modest gains in the AdaBoost model.

The 1:1 scatter plots in Figure 6 show that MS-1 yielded the best alignment between the predicted and observed SOC values, with data points distributed more symmetrically around the 1:1 line and a reduced bias observed in both the high- and low-SOC regions. In contrast, MS-3 showed marked underestimation at higher SOC levels and greater residual variance. The density-based heatmaps further support these findings: under MS-1, the predicted values were tightly clustered along the diagonal, indicating high agreement between the predictions and observations. This clustering decreased under MS-2 and was lowest under MS-3, where predictions deviated significantly, especially at high SOC levels.

In summary, the inclusion of VIs, particularly during periods of lush growth, substantially improved the predictive performance of the SOC models. MS-2, which incorporated VIs_bare, clearly proved more robust than the baseline MS-3 did, whereas MS-1, which used VIs_lush, emerged as the most scientifically grounded strategy for capturing cropland SOC variability. Among the three models, XGBoost made the most effective use of vegetation-related features, followed by the RF model, whereas AdaBoost was the least sensitive to VI input but still benefited from its inclusion.

3.2. Feature Analysis Based on the SHAP

Following the model performance evaluation, the XGBoost model was identified as the optimal model for determining the SOC content. To investigate the model’s internal decision logic, SHAP analysis was employed to interpret the contributions of individual features under MS-1, MS-2, and MS-3, as shown in Figure 7.

As shown in Figure 7, the SHAP summary plots under the three MSs demonstrate generally consistent rankings among most non-VIs. Key features, such as precipitation, DEM, SRD, and shortwave infrared bands (e.g., B7 and B5), consistently emerged as the top contributors across all the strategies. This consistency indicates the robustness and reliability of the model outputs under varying feature construction strategies. In contrast, differences were primarily observed in the ranking and presence of VIs, such as NDVI_lush and EVI_lush in MS-1 and their corresponding NDVI_bare and EVI_bare values under MS-2. These changes reflect the influence of temporal and contextual variations in vegetation information on SOC prediction, whereas the stable ranking of physical and topographic variables supports the validity of the adopted modelling strategies. Therefore, the relatively stable SHAP value distributions of non-VI features confirm that the ELMs are not overly sensitive to feature engineering choices, enhancing the interpretability and credibility of the inversion results.

Under MS-1 (Figure 7(a)(d)), VIs_lush and other vegetation indices were among the top-ranked predictors, demonstrating strong correlations with the observed ecological patterns. These variables showed broad and dispersed SHAP distributions—for example, NDVI_lush ranged from approximately –2 to +2, with both positive and negative contributions. This suggests a strong and nonlinear regulatory role in SOC predictions, which aligns with ecological processes such as photosynthetic carbon input and litter deposition [32,78]. EVI_lush highlights the multidimensional influence of vegetation structure and function during peak biomass periods. Notably, in the MS-1 strategy, NDVI_lush and EVI-lush contribute significantly to the explanation of SOC, while the SHAP value of GCI_lush is zero. We speculate that this is due mainly to its high correlation with NDVI_lush and EVI_lush. In XGBoost-based tree models, when multiple features are highly correlated, the model tends to prioritize more representative or informative features for node splitting, which results in the contribution of related features being diluted or replaced. Therefore, although GCI_lush itself may contain information, owing to its high overlap with NDVI_lush and EVI_lush, the XGBoost model prioritizes NDVI_lush and EVI_lush when the tree structure is constructed, resulting in the SHAP value of GCI_lush reaching zero. In contrast, MS-2 (Figure 7(b)(e)), which included VIs_bare, exhibited narrow and centralized SHAP distributions, mostly within the interval of –1 to +1. Their marginal contributions were limited, reflecting the lack of vegetation information and reduced spectral heterogeneity during bare-soil periods. This restricts the ability of the model to capture key ecological drivers of SOC, such as biomass input or microbial activity. The MS-3 strategy (Figure 7(c)(f)), which excluded VIs_lush and VIs_bare and relied on topographic and climatic variables (e.g., DEM, slope, and precipitation). The absence of phenological indicators reduced the ecological sensitivity of the model and impaired its ability to capture biologically driven SOC variation.

From an ecological perspective, MS-1 provides a more complete representation of vegetation–topography–carbon cycle coupling. It involves the integration of more rational VIs, spectral bands, and environmental attributes, enabling the model to depict the spatial heterogeneity of SOC more accurately. In contrast, MS-2, which is dominated by bare soil reflectance and static terrain factors, lacks the biological cues necessary to fully characterize SOC processes. Consequently, its explanatory power and ecological depth are significantly constrained.

In summary, SHAP-based analysis highlights the superior predictive capacity and ecological relevance of vegetation indices derived from the lush period. These features enhance model interpretability and support more accurate and ecologically meaningful SOC mapping by incorporating seasonal vegetation dynamics during periods of high productivity.

3.2. Cropland SOC Content

On the basis of the evaluation results from the test datasets, the XGBoost regression model was identified as the optimal algorithm for SOC content prediction. Consequently, it was employed under both modelling strategies, MS-1, MS-2 and MS-3, to conduct regional-scale inversion and perform statistical analysis of the predicted SOC content across the entire study area (Table 4).

Via the use of the optimal modelling combination—the XGBoost model under MS-1—the spatial distribution of the SOC content across croplands in Liaoning Province was mapped (Figure 8). The results reveal a distinct spatial gradient in the SOC concentration, reflecting the combined effects of topography, climate, land use, and vegetation productivity across the region.

As shown in Figure 8, the SOC content ranged from 2.19 g·kg⁻¹ to 33.86 g·kg⁻¹, with a mean value of 13.08 g·kg⁻¹, indicating moderate organic matter levels overall. The distribution pattern clearly displays spatial heterogeneity. Eastern Liaoning exhibits the highest SOC concentrations, particularly in the hilly and mountainous regions of Dandong, Benxi, and eastern Anshan. These areas are marked in red to orange in the map and are characterized by high precipitation, dense vegetation cover, and limited anthropogenic disturbance—all favourable conditions for organic matter accumulation. The SOC concentrations here often exceed 25 g·kg⁻¹, indicating high soil fertility. Central Liaoning, including parts of Shenyang, Liaoyang, and western Anshan, has moderate SOC levels (approximately 10–20 g·kg⁻¹), depicted in light green to light blue. This region includes major agricultural zones with a mix of irrigated cropland and hilly terrain, where SOC levels are influenced by both land use intensity and slope-driven erosion. In western Liaoning, particularly around Chaoyang and Fuxin, the lowest SOC concentrations, mostly less than 10 g·kg⁻¹, are shown in dark blue in the map. These areas are generally drier, more degraded, and experience frequent soil erosion, leading to limited organic matter accumulation and lower soil productivity. Southern coastal areas, such as Dalian and Jinzhou, exhibit moderate to low SOC levels, which is likely due to intensive land use, saline soils, and urban expansion. Overall, the spatial pattern aligns well with known climatic gradients, with higher SOC in humid and vegetated eastern regions and lower SOC in arid and erosion-prone western zones. The map clearly delineates agroecological zoning and can inform soil fertility management, land degradation monitoring, and targeted conservation practices.

These results also validate the effectiveness of MS-1, in which vegetation indices from the lush season and environmental covariates are integrated, in capturing the spatial variability of SOC across complex landscapes.

4. Discussion

This study provides a refined approach to estimate the SOC content in croplands using remote sensing and ELMs. Through model optimization, innovative feature selection, and comparative analysis, the proposed methodology demonstrates both theoretical significance and practical applicability. The following discussion elaborates on the methodological advances, model performance, and inherent limitations.

4.1. Innovation in Feature Selection and Modelling Strategy

The superior performance of MS-1 compared with MS-2 and MS-3 (Section 3.2) indicates that replacing bare-soil VIs with those derived from the crop-lush period is a key factor driving the increase in the SOC inversion accuracy. Previous studies have frequently relied on Vis, such as the NDVI and EVI, in SOC modelling, but most have relied on indices extracted from the bare-soil period [3,16,26,42,79,80,81]. However, such indices often neglect ecological dynamics during vegetation growth periods, which are more indicative of organic matter input. In this study, we proposed a novel strategy involving replacing bare-soil VIs (VIs_bare) with those from the crop-lush period (VIs_lush), guided by pedological theory on SOC formation. This approach yielded notable improvements in model performance. The vegetation density during the lush period is positively correlated with SOC, as denser vegetation generally contributes more plant litter and organic residues [82,83], thereby increasing the potential for the accumulation of SOC. VIs_lush, which reflect cumulative vegetation productivity during crop-lush periods, thus acts as an indirect proxy for SOC, integrating long-term effects of SOC on plant performance [84]. In contrast, VIs_bare primarily capture transient surface conditions and are more susceptible to non-SOC factors, such as soil moisture and roughness, which explains the lower accuracy of MS-2. Nonetheless, the reliability of VIs_lush is partially constrained by anthropogenic influences, such as fertilization and crop management [85]. Thus, while the use of VIs_lush enhances the predictive power, their effectiveness is context dependent and may not be universally applicable. Future research should assess this approach across different land use types to validate its robustness.

The SHAP dependence plots (Figure 9) further reveal the differences in contribution between bare-soil and crop-lush period VIs. NDVI_lush and EVI_lush exhibited clear and positive SHAP values, reflecting strong and consistent relationships with the SOC content. In contrast, their bare-soil counterparts showed weaker and more scattered relationships, indicating reduced predictive ability. Notably, GCI_lush exhibited SHAP values close to zero across the range of index values. Although the chlorophyll content, particularly as captured by GCI_lush, is mechanistically related to SOC through its link to photosynthetic capacity and biomass production, its predictive role in our models was minimal. This is likely due to the high correlation of GCI_lush with NDVI_lush and EVI_lush (Figure 10), resulting in information redundancy; consequently, the model prioritized NDVI_lush and EVI_lush for node splitting during training.

These findings underscore the importance of considering feature intercorrelations when modelling strategies are designed, as highly correlated variables may be downweighted despite their ecological relevance. Overall, our findings highlight that vegetation indices from the crop-lush period not only provide stronger and more interpretable signals for SOC prediction than their bare-soil counterparts do but also capture essential ecological processes underpinning SOC formation that static terrain or soil properties alone cannot represent.

4.2. Model Performance and Remote Sensing Data Considerations

Among the three ELMs tested, the XGBoost model consistently outperformed the RF and AdaBoost models following appropriate data preprocessing and hyperparameter tuning. Its iterative optimization mechanism effectively reduces prediction variance and handles structured datasets [68]. The performance hierarchy observed (XGBoost > RF > AdaBoost) agrees with findings from related studies [25,80], including recent SOC inversion research conducted in Liaoning Province within Northeast China, thereby reinforcing the reliability of the results of this study [44,79].

Methodologically, the dominance of the XGBoost model conforms with its ability to model nonlinear relationships and complex feature interactions, making it well suited for capturing the complex feedback loops among SOC, vegetation growth, and VIs_lush. This highlights the utility of ensemble learning in disentangling the multivariate drivers of SOC variability. However, model efficacy remains context-specific and is influenced by dataset characteristics, including class imbalance and noise levels. With respect to data sources, the accuracy of SOC inversion is also influenced by the spatial, temporal, and spectral resolution of remote sensing imagery. High-spatial-resolution platforms (e.g., Landsat and Sentinel-2) often suffer from low revisit frequencies, causing seamline mismatches in large-area mosaicking. Integrating multitemporal data can help mitigate these limitations and increase model stability.

Moreover, hyperspectral data—offered by several platforms, including Gaofen-5 and EnMAP—hold great promise for enhancing the accuracy of SOC prediction. Laboratory and field-based studies have validated the spectral sensitivity of SOC, and future efforts should prioritize the fusion of multispectral and hyperspectral observations for large-scale SOC content mapping [14,86].

4.3. Study Limitations

The proposed crop-lush period VI modelling strategy was evaluated primarily in cropland settings. However, in agricultural systems, SOC–vegetation relationships are often distorted by human interventions [52], such as fertilization and residue incorporation [55]. These factors can mask the natural feedback loop between vegetation and SOC. In contrast, forest ecosystems, which rely on natural processes, such as litterfall and root decay, may benefit more from crop-lush VIs in SOC inversion [32,33,87]. Thus, future work should include forest sampling to assess the transferability of this approach. In terms of the five principal soil-forming factors [88], incorporating temporal dynamics into SOC inversion poses a significant challenge. Accurate SOC content inversion requires the integration of both spatial heterogeneity and temporal dynamics. While topography and meteorological variables were incorporated in our approach, it did not account for temporally variable factors such as land use change and soil formation. Addressing these aspects could significantly enhance model reliability. Furthermore, the current model does not differentiate between paddy and upland fields—land use types with distinct carbon cycles. The incorporation of high-resolution land cover classification and time-sensitive predictors may reduce spatial uncertainty and improve prediction accuracy across diverse agroecosystems.

In this study, approximately 468 sample points were initially collected, and their uneven spatial distribution limited the generalizability of the model. To compensate, we applied kriging interpolation to augment the dataset to 2,799 samples. Although this increased spatial coverage, it also introduced interpolation-based uncertainties that may affect model fidelity. Future research should expand field sampling to establish a more robust training dataset.

5. Conclusions

In this study, we proposed an optimized modelling strategy for regional-scale SOC content estimation by integrating machine learning and remote sensing and recommended the replacement of bare-soil VIs with those from the crop-lush period. A cropland case study in Liaoning Province, Northeast China, demonstrated that distinguishing directly observed from indirectly influencing factors and prioritizing VIs_lush over VIs_bare significantly increased the accuracy of the SOC inversion. Among the three ELMs, the XGBoost model showed the highest sensitivity to VIs and achieved the best performance under MS-1 (R² = 0.84, RMSE = 2.22 g/kg, RPD = 2.49, and RPIQ = 3.25), followed by the RF model, whereas AdaBoost was less sensitive but still benefited from VIs_lush. SHAP analysis for MS-1, MS-2, and MS-3 confirmed the critical role of VIs_lush, revealing that vegetation conditions during the crop-lush period capture essential SOC variability beyond static terrain or soil features. The inversion results indicated an east–west SOC gradient, with higher values in the eastern mountainous and coastal areas and lower values in the central and western sandy regions (mean: 13.08 g/kg; range: 2.19–33.86 g/kg). This method supports long-term, large-scale SOC monitoring and provides valuable references for soil surveillance, crop planning, precision fertilization, sustainable land management, soil health assessment, and carbon cycling research.