Estimating Soil Organic Carbon Using VIS/NIR Spectroscopy with SVMR and SPA Methods

2.1. Study Areas

The study areas were located in Yixing (119°31′–120°03′E, 31°07′–31°37′N), Zhongxiang (112°15′–112°30′E, 31°15′–31°30′N) and Honghu (113°07′–114°05′E, 29°39′–30°02′N) (Figure 1). Yixing is located in the south of Jiangsu Province, China, with an annual temperature of 15.7 °C and a mean annual precipitation of 1177 mm [23]. Zhongxiang and Honghu are located in the middle of Hubei Province, China. Their annual temperatures are 15.0 °C and 16.6 °C with mean annual precipitations of 961 mm and 857 mm, respectively. The three areas are typical agricultural regions with various crops under cultivation.

2.2. Soil Samples

A total of 298 soil samples were collected from three study areas: Yixing (Figure 1b, 100 soil samples), Zhongxiang (Figure 1c, 100 soil samples) and Honghu (Figure 1d, 98 soil samples). The 298 soil samples were distributed over a wide range of land-use types, and most of the samples were taken from paddy field and rape field (Table 1). According to the Harmonized World Soil Database (version 1.2) [35], the dominant soil types in Yixing are anthrosols, luvisols, planosols and alisols for different crop cultivations; the soils collected from Zhongxiang mainly belong to anthrosols and leptosols, and those from Honghu belong to anthrosols and gleysols for planting crops (Table 2). For each soil sample, five soil subsamples taken within 1 m² were collected and thoroughly mixed to get a representative sample. At each sample site, about 1.0–1.5 kg of surface soils from a depth of 0–10 cm were collected after wiping off plant materials, plant residues, roots and stones. The collected soil samples were then kept in a sealed package for later spectral measurements and SOC content analyses in a laboratory.

2.3. Spectra Measurements and SOC Content Analyses

All the soil samples were air-dried at an indoor temperature for three days to standardize the moisture level, and small stones and plant residues were removed. The 298 soil samples were ground with an agate mortar and passed through a 20-mesh grid sieve (<2 mm). The reflectance spectra of all soil samples were measured in a laboratory using an ASD FieldSpec^®3 portable spectroradiometer (Analytical Spectral Devices, Inc., St, Boulder, CO, USA) with a wavelength range of 350–2500 nm, and the spectral readings were acquired in 1-nm increments over the wavelength range [36]). Each soil sample was put on a 10-cm diameter petri dish with a thickness of about 15 mm. Considering the optimal values of beam angle, lamp distance and sensor distance recommended by Zhou et al. [37], the spectra of all soil samples were measured in a dark room with a 50-W halogen lamp as the light source, which was positioned 0.3 m away from the soil sample, with a 15° zenith angle; while the optical probe was installed about 0.15 m above the soil sample. The correction with a standardized white Spectralon® panel with 100% reflectance (Labsphere, [38]) was made prior to the first scan and after every six samples. An average value of 10 spectral measurements for each sample was calculated as the final spectral reflectance. After spectral measurements, the SOC contents of all soil samples were determined using the Walkley and Black method [39].

2.4. Data Preparation

The SOC content measurements were combined with the spectra data of soils into a raw dataset. With principal component analysis (PCA) [40], the samples with both a large score distance and a large orthogonal distance were then detected as outliers and eliminated from the raw dataset. A detailed description of PCA can be found in Verboven and Hubert [40].

The reflectance spectra were reduced to 410–2450 nm in order to eliminate the noise at the edges of each spectrum [33]. The dataset of reduced spectra with SOC contents was then sorted based on the SOC content values from low to high, and the samples sorted at n and 2n orders (n ranging from 1 to 96) were selected as the calibration dataset, with 3n order samples as the validation dataset. Such a division ensures that the full range of SOC contents are represented in both the calibration and validation datasets [41].

2.5. Successive Projections Algorithm (SPA)

The successive projections algorithm (SPA) is a forward selection method, and its purpose is to select wavelengths whose information content is minimally redundant [26]. In the process of SPA, a projection operation in an orthogonal vector space is employed to select subsets of variables with a minimum collinearity. Variable selection by SPA is based on the principle that a new selected variable is the one among all the remaining variables having the maximum projection value on the orthogonal sub-space of the previous selected variable. The spectral data used for SPA wavelength selection are denoted by a matrix, X_{N × K}. Let M = min (N − 1, K) be the maximum number of variables to be selected. There are two sections to variable selection by SPA. The detailed SPA steps in the first section are described for a given initial wavelength, k(0), as follows.

• Step 1: Before the first iteration (n = 1), let x_j be the j-th column of the spectra matrix (X_{N × K}), j = 1,..., K.

• Step 2: Let S be the set of wavelengths that have not been selected yet. That is,

  $$\text{S}=\{\text{j},1\le \text{j}\le \text{K},\text{j}\notin \{\text{k}(0),\dots ,\text{k}(n-1)\left\}\right\}.$$

  (1)

• Step 3: Calculate the projection of x_j on the sub-space orthogonal to x_k(n−1) as:

  $$P{X}_j=X_j-\left(X_j^T{X}_{k(n-1)}\right)X_{k(n-1)}{\left(X_{k(n-1)}^T{X}_{k(n-1)}\right)}^{-1}$$

  (2)

where P is the projection operator.

• Step 4: Let k(n) = arg (max‖ Px_j‖, j∈S).

• Step 5: Let x_j = Px_j, j∈S.

• Step 6: Let n = n + 1. If n < M, go back to Step 1.

• End: The resulting wavelengths are {k(n); n = 0,…, M}.

The second section of SPA is to evaluate the candidate subsets of the variables selected in the first section. A total of M × K subsets of the variables are tested. The best variable subset is selected based on the smallest root mean squared error of prediction in the validation set of the MLR calibration.

In this study, the SPA was implemented using MATLAB version R2008b software (The Math Works, Natick, MA, USA) with a graphical user interface provided by Araújo et al. [43]. The graphical user interface for SPA is available at [44].

2.6. Model Development

Using the calibration dataset, the SPA-selected wavelength variables based on the eight pre-processed spectra were used to calibrate SVMR models, with a Gaussian radial basis function as the kernel function. The parameters, C and γ, of the kernel function were acquired by a grid-searching technique and a leave-one-out cross-validation procedure. The optimal values of these two parameters were selected when the minimal RMSE_CV was produced. The final SPA-SVMR model with an optimal pre-processing method was chosen.

In order to demonstrate the performance of the final SPA-SVMR model, a SVMR model without SPA operation and a SPA-PLSR model were developed with the same pre-processed spectra as the SPA-SVMR model. Then, these three calibration models were employed to estimate the SOC contents of the validation dataset. All models were implemented using PLS Toolbox 7.0 in MATLAB version R2008b software (The Math Works, Natick, MA, USA).

2.7. Model Evaluation

The coefficients of determination between the measured and estimated values for calibration (R²_C) and validation (R²_V) datasets, the error statistics root mean square error (RMSE) and the residual prediction deviation (RPD) [9] were used to evaluate the performance of the above models. The six level interpretations of RPD given by Viscarra Rossel et al. [45] were adopted as follows: RPD < 1.0 indicates very poor models/predictions, and their uses are not recommended; 1.0 < RPD < 1.4 indicates poor models/predictions, where only high and low values are distinguishable; 1.4 < RPD < 1.8 indicates fair models/predictions, which may be used for assessment and correlation; 1.8 < RPD < 2.0 indicates good models/predictions, where quantitative predictions are possible; 2.0 < RPD < 2.5 indicates very good, quantitative models/predictions; and RPD > 2.5 indicates excellent models/predictions.

3.1. Soil Organic Carbon Contents

According to the score diagnostic results of PCA analysis, a total of 11 samples with a large score distance and orthogonal distance were determined as outliers and were removed from the original 298 samples. Thus, 287 samples were left in the whole dataset, from which 192 samples were chosen as the calibration dataset and 95 samples as the validation dataset. The statistical descriptions of the SOC contents of the whole dataset, the calibration dataset and the validation dataset are presented in Table 3. Compared with the range of SOC content (0.79–30.73 g·kg⁻¹) for both the whole dataset and the calibration dataset, the validation dataset has a narrower range with 1.96–26.23 g·kg⁻¹, due to insufficient soil samples. However, the characteristic statistics of SOC contents in the calibration dataset and validation dataset were similar to those of the whole dataset; thus, the SOC contents of the calibration and validation dataset sufficiently represent those of the whole dataset.

3.2. Resample Intervals

As most wavelength information is redundant in hyperspectral data and the redundant information can decrease the efficiency of variable selection and modeling, proper resampling of hyperspectral data can be useful for rapid analysis and estimation of SOC contents. The resampling performance with 2−10 nm intervals was displayed in Figure 2. Considering the leave-one-out cross-validation results of SVMR and PLSR models calibrated using the calibration dataset, 9 nm was the best resample interval for both SVMR (R²_CV = 0.74, RMSE_CV = 2.78 g·kg⁻¹, RPD_CV = 1.95) and PLSR (R²_CV = 0.68, RMSE_CV = 3.12 g·kg⁻¹, RPD_CV = 1.74) models. After spectral resampling with 9 nm, only 227 variables were left for the subsequent variable selection.

3.3. Effects of Pre-Processing Methods on SPA and SPA-SVMR

The number of selected wavelengths (ranging from 9 to 186) by SPA and the distribution of these wavelengths varied strongly depending on the different pre-processing methods used in this study (Figure 3, Table 4). For example, SPA with SD selected the maximum number of wavelengths (186), and these wavelengths were located evenly throughout the entire VIS/NIR region, while SG selected 28 wavelengths, which were mainly located beyond the near-infrared region. Thus, these tests indicated that different pre-processing methods had a great influence on the SPA wavelength selection.

As compared with the SPA-SVMR model without pre-processing (Table 4), Log (1/R) + SG, SG, MC and MSC improved the model performance. SG achieved the best model performance (R²_V = 0.73, RMSE_V = 2.78 g·kg⁻¹ and RPD_V = 1.89), followed by Log (1/R) + SG (R²_V = 0.70, RMSE_V = 2.93 g·kg⁻¹ and RPD_V = 1.79). In contrast, the SPA-SVMR models with Log (1/R), FD, SD and SNV resulted in lower estimation accuracies, and Log (1/R) produced the worst performance (R²_V = 0.52, RMSE_V = 3.68 g·kg⁻¹ and RPD_V = 1.43). Comparisons of these tests demonstrate that different pre-processing methods make SPA-SVMR deliver distinctive performance results. Such results can be explained by the fact that several pre-processing methods (e.g., SG, Log (1/R) + SG, MC and MSC) removed noise and enhanced signals and, hence, were beneficial for the selection of feature bands by SPA. Thus, these pre-processing methods can improve the SVMR model accuracy.

3.4. Performance of SPA-SVMR, SVMR and SPA-PLSR Models

Table 5 shows the calibration and validation results of the SPA-SVMR, SVMR and SPA-PLSR models using the same SG-preprocessed dataset. The calibrations of three models achieved RMSE_V values ranging from 2.78 to 3.23 g·kg⁻¹. When compared with the SVMR model with 227 wavelengths (R²_V = 0.72, RMSE_V = 2.83 g·kg⁻¹ and RPD_V = 1.86), the SPA-SVMR model with only 28 wavelength variables achieved a similar estimation accuracy (R²_V = 0.73, RMSE_V = 2.78 g·kg⁻¹ and RPD_V = 1.89). The performance of the SPA-SVMR model was much better than that of the SPA-PLSR model, as R²_V and RPD_V were increased by 17.7% and 16.0%, respectively, while RMSE_V was decreased by 0.45 g·kg⁻¹. Therefore, the SPA-SVMR model was optimal for the estimation of SOC contents, since only 12.3% wavelengths were involved to achieve a similar or better estimation accuracy when compared with the SVMR and SPA-PLSR models.

4.1. Influence of Pre-Processing Methods on SPA and SPA-SVMR

The pre-processing methods used in this study showed considerably different influences on SPA when selecting wavelengths (Figure 3). One possible reason is that the collinearity between wavelengths varied in relation to the different pre-processed reflectance spectra in different degrees (Figure 4). As the main point of SPA is to obtain a subset of variables with minimal collinearity for SOC estimation, SPA with different pre-processing methods had strong differences in the number of selected variables and their distributions (Figure 3, Table 4). For example, after SD processing, collinearities among the wavelengths were considerably lower; therefore, SPA selected 186 variables out of 227 wavelengths. SNV generates similar collinearities for wavelengths as MSC (Figure 4f,h). Thus, the SPA-selected variables for SNV (20 variables) are similar to those of MSC (27 variables), except that there are a few different variables in the visible region (410−700 nm) (Figure 3). Hence, considering the collinearity between wavelengths of different pre-processed reflectance spectra, the variables selected with SPA vary in different degrees.

Based on the wavelengths selected by SPA with the eight pre-processing methods, eight SVMR models were correspondingly calibrated and validated (Table 4). SG can reduce the random noise and improve the quality of raw spectral data [42], which assisted SPA in selecting informative variables. Thus, the corresponding SPA-SVMR model achieved a better accuracy. Log (1/R) + SG also assisted SPA-SVMR to yield a good accuracy, but using Log (1/R) alone made SPA-SVMR deliver the worst performance among the eight pre-processing methods. Log (1/R) can reduce the nonlinearities, which is conflicted with the functionality of SVMR. Comparing with the SPA-SVMR without any pre-processing, SD and SNV obviously decreased the accuracy, because SD potentially generates noise [14] and SNV is sensitive to the noise in the spectrum [34]. Therefore, not all the pre-processing methods are appropriate to improve the accuracy of SPA-SVMR models. In this study, SG was an appropriate pre-processing method for the SPA-SVMR model. Moreover it was also employed and demonstrated to be useful in some other research [9,33,45–49]; for example, Vasques et al. [9] employed SG to improve the estimation accuracy for soil carbon.

4.2. Wavelength Selection

In this study, almost all of the 28 important wavelengths selected by SPA were located in the near-infrared region. Moreover, the eight pre-processing methods shared several wavelengths, particularly near 800, 1000, 1400 and 1900−2450 nm, which confirmed that the wavelengths in the near-infrared spectral region were important for SOC content estimations [4,16,22,24,30,50]. Ge et al. [13] estimated SOC contents with soil VIS/NIR spectra from different instruments/scanning environments. Their results also showed that most of the significant wavebands for SOC content estimation were located in the near-infrared region and shared by all the calibration models in their study. Such a result can be explained by the overtones and combinations of fundamental vibrations related to soil components; the stretching and bending of NH, OH and CH groups is concentrated in the near-infrared region [18,51]. Yang et al. [30] found that the visible range over 600−760 nm could be used to estimate SOC contents. However, this spectral region was not selected by SPA with SG in our study. A possible cause was the stronger influence of water and CH related bands, which overshadowed the influence of visible bands.

The selected wavelengths mainly located in the near-infrared region achieved the optimal estimation accuracy for SOC. Previously, some researchers reported that SOC estimations depend considerably on the spectral information over wavelengths greater than 1200 nm; this information might improve the accuracy of SOC content estimations [24,50,52]. Vohland and Emmerling [24] illustrated that the spectral features over 1900 nm selected by a genetic algorithm were of great significance for SOC content estimations. These results also agreed with the findings in this research, since 24 (85.7%) of the 28 spectral variables selected by SPA with SG pre-processing were located at wavelengths greater than 1200 nm. Specifically, the important wavelengths identified by SPA with SG for SOC coincided with those related to water (e.g., 1915, 1455, 1380, 1135, 940 nm) [16], organics (2381, 2147, 2084, 1967, 1931, 1517, 1454, 860 nm) [16,22], as well as clay minerals (2444, 2201, 2174 nm) [22], with a smaller number of wavelengths coinciding with those related to the iron oxides (e.g., goethite (920 nm), hematite (884 nm)) (Figure 3). These results were consistent with the findings of Viscarra Rossel et al. [16], who reported that most of the important wavelengths for estimating SOC contents were related to the absorptions of iron oxides, water, organics and minerals.

4.3. Comparison of SPA-SVMR, SVMR and SPA-PLSR Models

According to Table 5, when compared with the SVMR model, the SPA-SVMR model only employed 28 wavelength variables, but achieved a similar accuracy for the estimation of SOC contents. Such a result can be attributed to the fact that SPA chooses the informative wavelengths with the minimal collinearity [26] and removes many uninformative spectral variables, which simplifies the SPA-SVMR model and increases its robustness. This simplification of the SVMR model is also helpful in the interpretation of the SVMR model. SVMR, in contrast, may take much time for training samples when there are a large number of spectral variables and samples. In this case, SVMR may not obtain the best results, due to the redundancy in the input spectral data [53].

The SPA-SVMR model was superior to the SPA-PLSR model (Table 5). Our results confirmed the good performance of SVMR for estimating soil properties as compared to other multivariate calibration models [16,54]. The basic reason is that SVMR is more appropriate for dealing with potential nonlinearity from instrument variations (lamp aging and sensor sensitivity) than PLSR [8,16]. In addition, the SVMR model is less sensitive to noise and outliers than PLSR [18]. A similar comparison was also made by Liu and He [21] for determining the organic acids of plum vinegar. Their research also demonstrated that an SPA-SVMR model was better than an SPA-PLSR model. Liu et al. [28] also reported that an SPA-SVMR model outperformed an SPA-PLSR model and other models for the determination of the soluble solids of beer.

When compared with the previous studies in estimating SOC contents by using VIS/NIR spectra (R²_V = 0.57–0.91) [9,16,24,30,32,50,51,55–58], the predictive performance of the SPA-SVMR model (R²_V = 0.73, RMSE_V = 2.78 g·kg⁻¹ and RPD_V = 1.89) in this study was moderate according to the evaluation standard of Viscarra Rossel et al. [45]. Stevens et al. [54] used a large-scale EU soil survey of about 20,000 samples belonging to eight land-use types with a greatly wider range of SOC contents (0.0 to 586.8 g·kg⁻¹) to predict SOC contents and reported good prediction results, with R²_V values from 0.67 to 0.89 and RPD values ranging from 1.74 to 2.88. Brown et al. [50] achieved a R²_V value of 0.87 with a similarly wide range of SOC contents of 0.0–536.8 g·kg⁻¹ from a global scale. With 1011 soil samples and an SOC concentration range of 2.3–55.8 g·kg⁻¹, Shepherd and Walsh [54] obtained a R²_V value of 0.80 for SOC content estimation. Moreover, Saiano et al. [7] estimated the soil carbon contents of 89 soils from a small and homogeneous area, Pantelleria Island, and achieved a considerably higher accuracy with a cross-validation R² value of 0.951 and an RPD value of 4.49.

One basic factor restricting the further improvement in prediction performance of the SPA-SVMR is the range of SOC contents (0.79–30.73 g·kg⁻¹) and the number and representativeness of soil samples (287, Yixing, Zhongxiang and Honghu regions) in our study. Therefore, adequate samples over these regions with a wider range of SOC contents are required to further improve the prediction performance of the SPA-SVMR model. In addition, a range of soil types and land use types (Table 1, Table 2) may be another factor for the SOC content estimation with moderate accuracy in our study. Many studies have shown that the estimation models for soil properties varied to a large extent when considering soils diversity in various soil types and land use types [23,55,59]. That is because soils with different soil types can have considerably different characteristics, such as in compositions and spectra, due to the differences of the parent materials. Meanwhile, various land use types might further result in variations in soil properties, such as variation in observed SOC values. Thus, soil heterogeneity increases and model performance can be affected.

This study was carried out using the ground soil samples and the VIS/NIR spectra under laboratory conditions, and SPA-SVMR obtained good results when estimating SOC contents. These results laid the foundation for quickly estimating SOC contents from in situ field spectroscopy in future studies. However, for the natural soils in the field, the estimation of SOC contents could be affected by some factors, such as soil surface characteristics (e.g., particle size distribution, physical and biogenic crusts), mixed vegetation, atmospheric conditions and the varying moisture contents of natural soils, compared to laboratory-based dried soils. Recently, Minasny et al. [60] successfully utilized an external parameter orthogonalization algorithm to remove the effect of soil moisture from NIR spectra for the calibration of SOC contents and reported improved calibration and prediction of SOC contents. Wu et al. [61] identified a range of wavelengths in the NIR region in which the first derivative of the reflectance spectra seemed independent of the moisture content of the soil samples. They suggested to only use these selected wavelength intervals to obtain moisture-independent estimates of SOC under field conditions. Therefore, we will explore the proper pre-processing and algorithms for field soil samples and spectroscopy in depth to remove the effects of moisture and other environmental factors. The feasible SPA-SVMR model coupled with effective methods for the removal of natural environment effects may facilitate a good method for the rapid estimation of SOC contents from in situ field spectra in the future.