Inversion method of organic matter content of different types of soils in black soil area based on hyperspectral indices

2.1 Study area

The study area is shown in Figure 1. This area, situated in the black-soil region in and around Sihe Town, Yushu City, Jilin Province, experiences a temperate continental monsoon climate. The winter climate is cold, long, and snowy, whereas the summer rainfall is abundant. The average annual temperature is 5.3°C and the average annual precipitation is approximately 580 mm. The study area is situated in the hinterland of the Songliao Plain. The terrain within the area is relatively flat, with slightly undulating undulations and an average elevation of approximately 186 m. The region is rich in natural resources, including ample surface water resources, fertile soil, and a wide range of soils, mainly black soils, meadow soils, and albic soils, which are suitable for its growth, making it an important black soil planting area in China. Crops in the area are harvested once a year, with corn and rice being the main crops. Corn planting area accounts for >70% of the total cultivated area in the study area, which is a unique condition for the development of agricultural production.

Figure 1

Location map of the study area and photographs of three selected sampling sites (a) Jilin Province, (b) part of Yushu, (c) sampling area, and (d) sampling sites within the quadrangle (five-point method).

2.2 Sample collection and analysis

The distribution of agricultural croplands in the study area was investigated. Sampling points and routes were predetermined at regular intervals of 500 m in each direction, considering the available map data on soil and land use types. The samples were collected in mid-April 2022, and all collection sites were agricultural fields. Figure 1 shows the locations of the 518 sample plots. Each sample plot was 30 m × 30 m in size. Five soil samples were collected from each plot, containing three soil types: black soil, albic soil, and meadow soil. During the data collection period, the crops were harvested; there were no crops or weeds on the surface, and only a small amount of straw remained in some of the fields, which facilitated the collection of soil samples for analysis. In this study, a five-point sampling method was used to collect 500 g of soil at a surface depth of 5 cm [31]. The distance between the sampling point and other features (roads, buildings, etc.) was ensured to be >100 m to ensure that the image pixel at the sampling point was a “pure soil pixel,” meaning the pixel predominantly contained soil with minimal interference from other surface features, such as vegetation, roads, or water bodies [32]. Simultaneously, a handheld GPS locator was used to record the geographic coordinates of sampling points. After the soil samples were collected, they were placed in an oven for 48 h at 60°C to remove fine root impurities and ground [33]. The soil samples were subjected to indoor spectral measurements and organic matter content determination after being sieved through a 0.15 mm sieve [34]. In this experiment, the spectral reflectance of the soil samples was measured using an ASD FieldSpec 4 spectrometer with a wavelength range of 350–2,500 nm, spectral resolution of 3 nm, and wavelength accuracy of less than 0.5 nm [35]. To avoid measurement errors caused by the light source, all samples were measured in a dark room, and a standard whiteboard was used for calibration to reduce the noise and maintain a consistent distance of 3 cm between the probe and the surface of the sample. The measurement was repeated ten times for each soil sample, with an interval of 1 s between each measurement. The final results of the ten measurements were averaged to obtain the measured spectral data on the ground. The SOMC of the soil samples was determined using an external heating method with potassium dichromate-concentrated sulfuric acid. This method is accurate and suitable for the analysis of a large number of samples because of the non-interfering effects of carbonates in the soil [36].

2.3 Data preprocessing

The ZY1-02D hyperspectral image data used in this study were provided by the China Center for Resources Satellite Data and Application and contained hyperspectral data in 166 spectral bands with a spectral coverage of 400–2,500 nm. For this study, the satellite image acquired on 26 April 2022 during a cloud-free observation window was selected to ensure synchronization with soil sample collection. During the period from sample collection to imaging, there was no agricultural production activity in the study area, and the change in SOMC under natural conditions was extremely weak, which did not affect the sample collection. The streaking phenomenon is evident in the ZY1-02D hyperspectral short-wave infrared (SWIR) data; therefore, the global debanding method is used to repair the streaks. For the effect of water vapor in the atmosphere, the satellite hyperspectral images are removed from the water vapor absorption bands at 1,350–1,450 nm and 1,800–1,900 nm wavelength range of the band that is seriously affected by water vapor, as well as the edge spectral channels containing a large amount of noise data before 400 nm and after 2,470 nm wavelengths [37]. A total of 21 spectral bands were cumulatively eliminated and the processed hyperspectral image contained 145 bands. The processed hyperspectral image contained 145 bands. Radiometric calibration was performed according to the satellite calibration file and the FLAASH model was used to correct the atmosphere of the radiometrically calibrated ZY1-02D hyperspectral data. To accurately identify the study area, the SVM classification method was used to supervise the classification of non-agricultural land parcels, and cloud shadows, roads, and residential areas in the images were eliminated (Figure 2).

Figure 2

SVM method for supervised classification extracted bare soil image element results.

2.4 Selection and construction of spectral index

Difference indices (DI), normalized indices (NDI), and ratio indices (RI) are commonly used data processing and analysis methods that are widely applied to extract relevant surface information [38]. The DI characterizes the spectral difference by calculating the difference in reflectance between two spectral bands; the NDI is used to normalize the difference between two bands; and the RI characterizes the relative difference between bands by calculating the ratio of the reflectance of two bands. These indices help highlight the spectral differences of specific features and can effectively eliminate the differences caused by light, sensors, atmosphere, and other factors, making the changes in feature characteristics more obvious, thus enhancing the interpretation and analysis of remote sensing data [39]. The calculation method of DI can be used to highlight the features of different features, RI can strengthen the reflectance differences of features in certain bands, and NDI can reduce scale differences. The formulae are as follows:

Most currently published spectral indices are derived from the physical characteristic laws of spectral reflectance, including the vibrations and combinations of reactive groups such as C–H groups, C–N groups, C═O groups, and –OH groups, and electromagnetic activity resulting from electron jumps occurring in the VN-SWIR spectral interval [40]. These spectral indices consider the interactions between electromagnetic radiation, atmosphere, vegetation cover, and soil background [41]. These have been gradually improved and developed through continuous improvements in mathematical simulation statistics. For example, Soil Adjusted Total Vegetation Index (SATVI) detects the reflectance of vigorous and senescent vegetation through the ratio of red and two SWIR radiations. In combination with the soil correction factor, this can reflect the nature of the soil in densely vegetated areas [42]. Brightness Index (BI) explains the relationship between soil brightness and remote sensing data through the mathematical operation of reflectance in the two visible bands and is more sensitive to the texture of the soil [43]. However, the relationship between spectral indices and organic matter content in bare soil areas has not yet been studied in detail. To investigate the potential of this remote sensing index in predicting the organic matter content of bare soil areas, this study was designed to develop two types of spectral indices: the vegetation index and the luminance correlation index. These indices were based on six types of spectral indices, using the waveband of ZY1-02D. The spectral index calculation formula and the corresponding waveband channels are presented in Table 1.

The Pearson correlation coefficient (PCC) is a statistical measure that quantifies the strength of the linear relationship between two variables, with values ranging from −1 to 1 [48]. For two-band indices, calculating the correlation coefficients between the spectral indices of the sample set and SOMC is the most commonly used and effective method for extracting the characteristic spectral indices. In this study, two-band indices were calculated based on ZY1-02D hyperspectral remote sensing images. The spectral indices DI, NDI, and RI of any two bands R _i and R _j were calculated by traversing the full-band soil spectral data from 400 to 2,500 nm. The combination of TVI and BI was restricted within the band channel range. PCCs between various indices and SOM were calculated, and a correlation coefficient matrix was constructed to analyze the characteristic spectral indices of SOM. Most existing studies on SOM estimation using spectral indices use two-band spectral indices for environmental modeling and attribute quantification. However, only a limited number of studies have explored the potential of three-band spectral indices to estimate soil attributes. For the three-band indices, the spatial scatter plots of the correlation coefficients between the three-band index and SOM were generated by calculating the combinations of spectral indices within the three spectral band channels [α₁,α_n], [β₁,β_n], and [γ₁,γ_n], which correspond to any three bands R_α, R_β, and R_γ. The methodological workflow of this study is shown in Figure 3.

Figure 3

Flowchart of research methods.

2.5 Model establishment and accuracy evaluation

In this study, the XGBoost algorithm was used to build the SOM prediction model. XGBoost is an efficient and scalable gradient boosting algorithm [49]. It is based on the concept of boosting trees, which integrates multiple weak learners (typically decision trees) to create a strong learner, thereby enhancing the model’s predictive performance. Each new tree is built by iteratively reducing prediction errors, allowing the model to better correct the mistakes of the previous models. To prevent overfitting, XGBoost incorporates regularization terms (L1 and L2 regularization) during the construction of each tree, effectively controlling the model’s complexity. After the model creates the booster tree, it is relatively straightforward to obtain the importance score for each attribute [50]. In general, the importance score measures the value of a feature during the construction of a boosting decision tree in the model. The more an attribute is used in the construction of the decision tree in the model, the higher its relative importance. The importance was calculated in a single decision tree based on the amount of improvement in the performance metrics by each attribute split point. Nodes are responsible for evaluating and recording the frequency with which an attribute carries more weight to enhance the split-point performance metric. The importance of an attribute increases as it is selected more frequently by decision trees. The evaluation metrics selected for this study were the coefficient of determination (R ²), normalized root mean square error (NRMSE), and Residual Predictive Deviation (RPD). The expressions for the above evaluation metrics are as follows:

where n denotes the number of samples in the dataset, y i is the true value, y ˆ i is the predicted value, and y ¯ i is the average of the true values.

3.1 Soil reflectance spectral characteristics

To analyze the spectral reflectance characteristics of different soils in the study area, the ground spectra of 518 soil samples collected were statistically analyzed by grouping them according to soil type. Subsequently, the measured indoor spectra and satellite spectra were compared separately. Figure 4(a) shows the average spectral reflectance of three soils with the same organic matter content. Overall, the spectral reflectance of meadow soil is higher, while the reflectances of black soil and albic soil are similar in the visible region; however, after 1,000 nm, the reflectance of black soil exceeds that of albic soil. In the near-infrared region, there are significant differences in the spectral reflectance of the three soils. The average spectral reflectances of the black soil, albic soil, and meadow soil at different SOMCs are represented by b, c, and d in Figure 4. In the wavelength range of 400–2,500 nm, the overall soil reflectance exhibited a decreasing trend with increasing SOMC. It is noteworthy that the differences in reflectance among the three groups become more pronounced in the NIR region, especially beyond 1,000 nm, where the reflectance of meadow soils is more affected by SOMC. A comprehensive comparison of the three sets of spectral data revealed that although the ZY1-02D hyperspectral data were affected by conditions such as soil roughness and soil moisture, the change rule with the indoor measured spectra was generally aligned with that of the pixel spectral reflectance in a specific wavelength range, with the variation in SOMC being more significant.

Figure 4

Comparison of average spectral curves of three soils under the same SOMC and different SOMC conditions.

To elucidate the relationship between the satellite spectra and indoor measured spectra, the correlation coefficients of the two spectra of the three soil-type samples with SOM were calculated. As shown in Figure 5, the spectral correlation coefficient curves of the three soil types indicate that SOMC primarily affects the VNIR spectral range, with the overall correlation ranked from high to low as black, meadow, and albic soils. However, in the range of 460–1,000 nm, the correlation of albic soil began to surpass that of meadow soil. The correlation coefficients of the three soil types in the visible region exceeded 0.4, indicating a strong correlation. This demonstrated that the spectrally corrected image element spectra retained the majority of soil spectral features and could be effectively used for SOM estimation.

Figure 5

Correlation between reflectance of two spectra (measured spectrum and ZY1-02D spectrum) and SOMC.

3.2 Spectral indices

To accurately extract the wavelength combinations with the highest correlation between SOM and dual-band indices, this study was designed to compute five spectral indices (DI, NDI, RI, TVI, and BI) for the combinations of two bands in the range of 450–2,500 nm, and to analyze the correlation between them and SOM. Figure 6(a)–(e) shows the five dual-band index and SOM correlation heat maps. The results indicated that the two strongest correlations of the band were primarily found in the visible and near-infrared wavelengths and were predominant at 1,000 nm or below. The high correlation with SOM in DI is the combination of the 400–580 nm and 580–860 nm bands, NDI is the combination of the 490–540 nm and 550–590 nm bands, and RI is the 490–540 nm and 550–580 nm. Although the correlation between TVI and SOM was generally low, the index still provided information that could reflect the nature of SOM in the spectral combination of 700–800 nm and 600–700 nm. To make full use of bands rich in characteristic information in the images, this study selected spectral indices of the band combinations with correlation coefficients >0.5 in the satellite spectral combinations for analysis. In total, 56 DIs, 40 NDIs, 36 RIs, 33 TVIs, and 15 BIs were screened.

Figure 6

Heat map of correlation between DI, NDI, RI, and SOM. (a) DI, (b) NDI, (c) RI, (d) TVI, and (e) BI.

A scatter plot of the correlation between the four types of three-band indices and the organic matter content is shown in Figure 7. The spectral band combinations with higher correlations between SATVI and SOM were mainly concentrated in the bands of 630–690 nm, 1,560–1,700 nm, and 2,080–2,300 nm, and the band regions with higher correlation of EVI were 470–510 nm, 600–640 nm, and 840–1,020 nm. Owing to the broad spectral bands of SATVI and EVI, there were up to 2,000 spectral band combinations. To minimize information redundancy, only the band combinations with correlation coefficients >0.5 were selected, resulting in 24 and 27 band combinations for SATVI and EVI, respectively. The spectral combinations of The Second Brightness Index (BI2) and Redness Index (RI) were generally moderately correlated with organic matter content, with a few combinations having weak correlations, and a total of 24 and 18 band combinations with correlations >0.5 were selected.

Figure 7

Degree of correlation between four three-band indices and SOMC.

3.3 Model construction for split-sample strategy

Considering the information redundancy of hyperspectral indices, XGBoost ranked the importance of the 273 spectral index features in the modeling. Feature variables were iteratively selected in descending order of importance, with optimal variable selection determined using K-fold cross-validation to obtain the model with the lowest NRMSE. Figure 8 shows the model NRMSE for all types of soil samples when cross-validation was used to select the optimal features. The model performance gradually improved with an increase in the spectral index. However, after the index increased beyond the optimal spectral index, the model performance tended to become weakened. The best inversion model accuracy was achieved by setting the retained feature variables to the first 27 in the sample groups of all soil types and meadow soils, with NRMSEs of 10.07 and 8.77%, respectively. The highest model accuracy was achieved by selecting the first 23 and 25 spectral indices for the black soil group and albic soil group, resulting in NRMSEs of 8.17 and 9.42%, respectively. Table S1 lists the optimal set of variables for the four soil types.

Figure 8

Iterative analysis of the number of features and XGBoost model performance for different sample groups.

3.4 Spatial estimation and mapping of SOMC

In this study, the optimal spectral indices for each group were used as independent variables in the model. SOMC was used as the dependent variable, the number of samples was randomly divided into training and validation groups at a 3:1 ratio, and local regression estimation was performed on different soil sample types using the XGBoost model. Table 2 presents the regression accuracies of the three SOM estimation models. The estimation model for black soil used 23 spectral indices as dependent variables, and the model’s coefficient of determination was 0.845. The estimation model for albic soil used 25 spectral indices as dependent variables, and the coefficient of determination of the model was 0.744. The estimation model for meadow soil used 27 spectral indices as the independent variables, and the model’s coefficient of determination was 0.788. This indicates that the estimation models under the split-sample strategy exhibited better performance.

The study area was divided into three soil subzones based on the soil type. The content was estimated and the distribution of SOMC was mapped using the XGBoost model. This was based on the characteristic spectral indices and extracted bare soil image. Figure 9 illustrates the spatial distribution of SOMC within the study area, which ranged from 0.4 to 2.89%, with an average SOMC of 1.24%. The overall SOMC level was predominantly 1.5–2%, accounting for 43.39% of the total study area, with the highest 2–3% range accounting for 15.73% of the study area. The spatial distribution of SOMC was generally high in the east-central region, while areas with lower values were primarily dispersed in the southern region. The high-value areas are roughly banded due to the cultivation of paddy fields in the eastern region, where long-term application of nitrogen fertilizers and irrigation results in a high humus content.

Figure 9

Spatial distribution of predicted SOMC under split-sample regression conditions.

4.1 Accuracy analysis of SOM estimation model based on hyperspectral index under split-sample strategy

To further analyze the superiority of the split-sample strategy, the R² and NRMSE values of the SOM estimation models were used to assess predictive accuracy across soil types, with black soil achieving the highest R² of 0.845 and the lowest NRMSE of 8.17%. Figure 10 shows the error scatter plots of the split-sample and full-sample regression models categorized according to the three soil types. The SOMC estimation model accuracies were in the following descending order: black soil, meadow soil, albic soil, and full-sample groups. The estimation model of the full-sample regression had 27 spectral indices as independent variables and the coefficient of determination of the model was 0.560. Among the regression models for the three soils, the R² of the SOMC estimation model was 0.845, 0.744, and 0.786 for the black, albic, and meadow soils, respectively. It was observed that the spectral index model of the full-sample group was limited by the complexity of the soil type and limited accuracy improvement. These results indicate that the models using the split-sample strategy demonstrated superior predictive ability and a better fit to the data. Among them, black soil had the best accuracy for the SOMC estimation model, which may occur for two reasons. First, it is caused by different physical properties of the soil structure. The structure of black soil is looser and mainly arranged in agglomerates or blocks, which helps to maintain aeration and drainage of the soil [51]. In contrast, the soil structures of the albic soil and meadow soil were mainly dominated by blocks or columns, leading to relatively poor soil drainage [52]. Particularly in albic soil, the albic layer has a high hardness and water capacity, inevitably causing water stagnation in the upper layer. Considering the varying intrinsic soil storage water content, this difference may result in variations in the accuracy of estimation among different datasets. Second, it may be caused by different chemical components of different soil types. Compared with albic soils, black soils and meadow soils have rich humus layers and higher SOMC accumulation, which increases the variability of the spectral signals. Lu et al. [53] found that the accuracy of the spectral prediction model for SOM was not satisfactory when considering traditional soil types. This is attributed to the presence of sand ginger black and yellow-brown soils, both of which contain high levels of clay minerals, particularly montmorillonite. The SOMC values in these two soil types were relatively similar, leading to spectral curves that exhibited similar characteristics. As a result, the regression accuracy for SOMC was low. In contrast, the clay minerals in the black and meadow soils in this study exhibited significant differences, which can enhance the acquisition of valuable information, thereby improving the accuracy of the model.

Figure 10

Comparison of accuracy between split-sample regression and full-sample regression based on three soil type classifications.

Therefore, by analyzing the spectral differences of different soils and then constructing a spectral index applicable to the assessment of different soil properties using methods such as correlation analysis and feature selection, it helps in soil quality evaluation. However, although the hyperspectral index developed in this study enhanced the differences in the independent variables and the accuracy of the estimation model, it still excluded a significant amount of information from the bands, which may have resulted in the omission of information and limited the accuracy improvement of the estimated model. To maximize the use of valuable information within the bands, the derivative method can be applied to hyperspectral data to highlight the distinctions among spectral information [8]. This approach can improve the estimation accuracy of the surface parameters while effectively capturing the fine details of the spectral curves.

4.2 Analysis of the contribution of eigenspectral index variables in regression models

Combined with the analysis of the spectral reflectance of various soil types, it was found that different soils exhibit different response characteristics in satellite spectral bands, as well as in different hyperspectral indices. To explore the degree of response of various soil types to these spectral indices, the importance scores of the spectral index features in the estimation models for the different soil types were assessed, as shown in Figure 11. This was performed by calculating and ranking each attribute in the dataset, as shown in Figure 11. This result was obtained from the XGBoost model by weighting and summing the results of each attribute in all boosted trees and then averaging them. This process highlights the differences in the spectral indices of each group within the sub-sample while also reflecting the varying degrees of contribution of different spectral indices in the model. In both the full-sample and split-sample regressions, the DI index, with R₄₉₉ and R₅₇₆ as the dependent variables, exhibited a significant influence on the prediction of the four model groups. The proportions of DI and NDI in the spectral parameters of the black soil estimation model reached 39%. Moreover, the proportions of DI and NDI were higher in the albic and meadow soil estimation models. In contrast, ratio indices (RI), although not a large proportion of the model and does not have the ability to significantly improve its predictive ability, may play a role in stabilizing the model when participating in the model construction process together with other spectral indices.

Figure 11

Characteristic scores of spectral indices involved in modeling under different soil types.

When distinguishing between the three soil types for analyzing the F-score, an increasing trend in the importance of vegetation indices in the model was observed. In particular, the EVI, which consists of the R₄₅₆, R₆₇₉, and R₉₃₇ bands, exhibited a higher degree of characteristic response to SOM across the three soil backgrounds. The average NRMSE of the model decreased by 0.2% at the second iteration in the EVI for the black soil group. Additionally, the average NRMSE decreased by 0.3 and 0.1% at the third and twenty-first iterations in EVI and SATVI for the albic soil group, respectively. These results indicate that the vegetation indices, such as EVI and SATVI, exhibit varying degrees of responsiveness to organic matter across different soil types, aiding in SOMC estimation in areas with complex soil types. This is consistent with He et al. [24], who found that the near-infrared band and SATVI are crucial for SOC prediction. This finding supports the notion that the vegetation index responds differently to organic matter in various soil types during the process of soil type differentiation, highlighting its potential to enhance the accuracy of estimating SOMC in bare soil areas.

4.3 Uncertainty analysis of modeling using hyperspectral indices

As demonstrated by the accuracy of the SOM estimation model, spectral indices serve as reliable inputs for the construction of estimation models. Specifically, spectral indices such as DI and NDI played a pivotal role in enhancing both model stability and prediction accuracy. They effectively captured critical variations in soil reflectance at key wavelengths, particularly at 499, 516, 559, and 576 nm, which are highly sensitive to organic matter content. This ability to comprehensively capture spectral responses related to SOM significantly contributes to the overall effectiveness of the estimation models.

In the optimal variable set for the black soil group, the combination of bands R₅₇₆ and R₅₈₅ was used most frequently, each accounting for 14%. For the albic soil and meadow soil groups, R₅₇₆ had the highest utilization rate. The correlation coefficients of these bands in previous correlation analyses exceeded 0.48. Therefore, it can be seen that the band combinations of spectral indices in various prediction models were primarily concentrated in the 500–600 nm range, which is consistent with the broad absorption bands identified by Angelopoulou et al. [54] in their spectral study of soil properties. Different types of soil contain different forms of humus with diverse chemical compositions. Weak overtones and combinations of these vibrations are present in NIR-SWIR (700–2,500 nm), owing to the stretching and bending of the N–H, O–H, and C–H bonds in the soil components. From the results of the iterative analysis, it is clear that vegetation indices such as EVI and SATVI have high feature scores in the regression model for distinguishing soil types. They help provide more comprehensive information about soil properties.

However, there are still some uncertainties in the estimation of organic matter content using spectral indices. Factors such as small quantities of straw, moisture, and atmospheric conditions can influence the accuracy of organic matter content estimation. Many of these variables change over time, which significantly exacerbates the uncertainty problem. These include temperature, precipitation, soil moisture content, and the physical canopy structure. Therefore, finding and synergizing important environmental variables with the specific conditions of the study area, and developing a bare soil spectral index that can correct for atmospheric and vegetation effects to optimize the estimation model is one of the important directions to improve the accuracy of SOM inversion tools in the future.