Assessment of Soil Characteristics Using a Three-Band Agricultural Digital Camera

The research was conducted within two arable fields located in Pokrzywno (Wielkopolskie Voivodeship, Poznań Poviat). This region has a temperate transitional climate characterised by a small number of frosty days and low rainfall. The average annual temperature is 8.5°C and the annual rainfall is about 500–550 mm (WIOS 2013). It is an area with unfavourable water balance, exposed to periodic droughts. Soils classified as Luvisols and Phaeozems, according to the IUSS Working Group WRB (2015), dominate this study area (Fig. 1).

Fig. 1

A total of 151 samples were collected from both research fields. They were tested in the soil science laboratory of the Adam Mickiewicz University in Poznań. All samples were prepared for testing by drying, grinding in a ceramic mortar and sieving through a 2-mm mesh sieve. The soil texture was determined by the hydrometer method according to the standard PN-R-04032 (Polish Committee for Standardisation PKN 1998). SOC was determined using oxidation by K₂Cr₂O₇ with H₂SO₄ for 30 min on the digestion block at 150°C and titration of oxidant residues by FeSO₄ (Nelson, Sommers 1996). Total nitrogen was modified using the Kjeldahl method (International Standard ISO 12261 1995). The soil pH was measured in 1:1 soil solution ratio in water and 1M KCL (PN-ISO 10390 1997). The form of nutrients available to plants (K, Mg, Ca, Zn, Cu, Pb, Cd, Mn and Fe) was determined by the modified Mehlich 3 method (Mehlich 1984). CEC was determined by successive barium and magnesium chloride solution extraction and flame atomic absorption spectroscopy (International Standard ISO 11260 1994). Calcium carbonate content was determined twice using Scheibler volumetric (International Standard ISO 10693 2002) and titration methods (FAO 2021). As a result of the analyses, the soil particle size composition was determined, and the amount of organic carbon, nitrogen, the ratio of carbon content to nitrogen content, soil reaction, percentage of calcium carbonate, total CEC and the content of the elements potassium, magnesium, calcium, zinc, copper, lead, cadmium, manganese, iron and phosphorus was also determined.

Multispectral data of soil samples were taken in the laboratory environment using the ADC by Tetracam. The specific design with optical Bayer filter mask in the complementary metal–oxide–semiconductor (CMOS) sensor allows to obtain three images at a resolution of 2048 × 1536 px (3.2 Mpx) (Swain et al. 2010). The images correspond to three Landsat Thematic Mapper 2, 3 and 4 spectral bands: green (520–600 nm), red (630–690 nm) and near-infrared (760–920 nm) (Lan et al. 2010) and the estimated ground pixel resolution is 0.000707 m px⁻¹ (Swain et al. 2010).

For the purposes of this study, photographs of 151 soil samples were taken under laboratory conditions. The ADC was placed on a tripod at a height of 70 cm above the test object and at an angle of 90°. In addition, a 400-W halogen lamp was used to illuminate the surface of the soil, which was set at a distance of 80 cm and at an angle of 45°. Then the images were processed into tiff format in the Pixel Wrench 2 program dedicated by the manufacturer. The next step was to transform the original digital numbers to reflectance using the TNT Mips software.

The evolution of remote sensing techniques caused the development of methods of evaluating remote measurements, processing and extracting as much information as possible from the collected data. Attempts to interpret the reflectance data from different available ranges of electromagnetic radiation have led to the development of a large number of indicators and their derivatives. A large group of spectral indices relate to the vegetation and soil substrate, which are calculated as the ratio of reflectance in two or more spectral channels of the selected device, sometimes with additional parameters (Bannari et al. 1995).

In this study, vegetation indices such as normalized difference vegetation index (NDVI), radar vegetation index (RVI, Bannari et al. 1995, Martínez M. 2017) and infrared percentage vegetation index (IPVI, Gunathilaka 2021) were used. All of them presented ratios between reflected radiation in red and infrared ranges. In addition, the indexes with green spectra such as IPVI (GNDVI, Candiago et al. 2015), green-red vegetation index (GRVI, Motohka et al. 2010) and modified GNDVI normal (Crippen 1990) were added to the dataset. Finally, three variants of the soil-adjusted vegetation index (SAVI, Bannari et al. 1995) were calculated, considering the soil parameter in three values of 0.25, 0.5 and 0.75. Summarises the spectral indices used in the spectral data processing (Table 1).

Raw spectra are subject to fluctuations and noise disturbance. For this reason, methods of standardising spectra data are often used. They consist in reducing the undesirable effects in the set of spectral measurements (Gholizadeh et al. 2015). In this study, we used methods of multiplicative scatter correction (MSC), standard normal variate (SNV), conversion of reflectance data into absorbance and scaling with minimum and maximum values. One of the most popular methods of data standardisation is MSC (Rinnan et al. 2009). This method relies on adjusting to each spectral measurement an ideal reference spectrum estimated based on additive and multiplicative correction factors (Rinnan et al. 2009). Another frequently used method is the SNV, which consists of common centring and scaling by subtracting the mean values and normalising with the standard deviation for each reflection spectrum (Vestergaard et al. 2021). MSC and SNV were introduced by the Prospectr 0.2.4 package implemented in the R-4.1.3 software for Windows. When working with various types of data, it can be noticed that occasionally, some measured results may significantly differ from others and thus disturb the work of the computational model (Gholizadeh et al. 2015). To eliminate this effect, the min–max scaling can be used (available in the R software in the Caret package). This type of data normalisation is based on scaling all data so that they fall in the range from zero to one. This reduces the value of the standard deviation and also the effect of outliers in the dataset.

Additionally, the spectral data were converted to the form of absorbance according to the following formula (Wenjun et al. 2014): ABS=−logR, ABS = - \log R, where R is the reflectance in a given spectral channel.

Regression models are used to establish the relationship between variables y (dependent, explained variable) and x (independent, explanatory variable). Such analysis makes it possible to explain how the value of the explained variable developed under the influence of the explanatory variable. In the case of a greater number of variables, it is useful to use statistical programs that calculate relationships between the variables by selecting the best variables to determine the relationship within the given parameter. Regression models are often used to estimate the content of soil components. Among them are Cubist (Peng et al. 2015) and PLS (Vestergaard et al. 2021). Cubist, which is modelled on the M5 Quinlana model (Quinlan Basser 1992), is a tool that allows to create a decision tree based on a given number of rules. Each rule represents a linear regression model to which the given variables are fitted. If the value of the variable matches the rule, its predicted value is calculated (Minasny, McBratney 2008). The final regression model is simplified to reduce the absolute error value (Yi Peng et al. 2015). Partial least squares (PLS), which was created by H. Wold in 1966, is a linear regression method recommended in the case of a large number of explanatory variables and when there is a high probability that they can be correlated with each other. The reduced set of x variables is used to create the regression model (Wold et al. 2001).

Variable importance in the projection (VIP) is useful for determining which predictor variables are best explained by explanatory variables. VIP determines the variables and the extent to which they contribute to the construction of a given regression model (Chong, Jun 2005, Xu et al. 2021).

VIP values can be obtained through dedicated software.

There are many ways to determine how well an outcome estimate is guaranteed by a given regression model. One of them is to calculate the coefficient of determination R². This measure shows how well the variance of the predicted spectral values coincides with the values measured in the laboratory. R² is calculated according to the following formula (Ng et al. 2022): R2=1−∑in(xi−yi)2∑in(xi−x¯i)2, {\rm{R}}2 = 1 - {{\sum\nolimits_i^n {{{\left( {{x_i} - {y_i}} \right)}^2}} } \over {\sum\nolimits_i^n {{{\left( {{x_i} - {{\bar x}_i}} \right)}^2}} }}, where x_i i y_i are the values observed and predicted under i, and is the mean value of the observed (Ng et al. 2022).

Another factor is the root mean squared error (RMSE), which informs about the difference between the values estimated in the model. RMSE takes values equal to or greater than zero, with zero being a statistically perfect match of estimated values to those observed (Peng et al. 2015): RMSE=∑in(xi−yi)2n, RMSE = \sqrt {{{\sum\nolimits_i^n {{{\left( {{x_i} - {y_i}} \right)}^2}} } \over n}} , where x_i i y_i are the values observed and predicted under i, and n is the number of samples (Yi Peng et al. 2015). Rel-RMSE – the derivative of RMSE is the relative mean square error calculated as the ratio of the RMSE value to the mean value of a given variable (Yiping Peng et al. 2019). Regression point displacement (RPD) is another value that is used to evaluate the quality of the model. It is the ratio between the standard deviation of the variable and the RMSE (Vestergaard et al. 2021). At present, the RPIQ measure is considered more frequently when assessing the model estimation. It considers only the interquartile range of the results. It is calculated as the quotient of the values measured in the interquartile space and the RMSE (Vestergaard et al. 2021).

Saeys et al. (2005) proposed to establish the criteria for the classification of the model according to the following values of R² and RPD: if R² is <0.5 or RPD is <1.5, this indicates a poor model estimation and means that the calculated values cannot be used. If R² is in the range of 0.5–0.65 and RPD is between 1.5 and 2, then it is possible to distinguish between high and low values. R² from 0.66 to 0.81 or RPD from 2 to 2.5 shows that the model enables approximate quantitative predictions. When the value of the coefficient of determination is in the range of 0.82–0.9 or RPD is between 2.5 and ≥3, then the model can be assumed to be good. Finally, if R² is >0.91 combined with value RPD >3, the model is perfect. As a rule, a good prediction model should have the highest R² i RPD (or RPIQ) values and as little RMSE as possible (or rel-RMSE) (Wenjun et al. 2014).

The results of laboratory analyses are presented in Table 2. The summary includes statistics on mean values, median, maximum and minimum values, and standard deviation. As shown in Table 2, mean values for all data are in the range of 0.12–1504.29 mg ∙ kg⁻¹. Calcium has the most varied values. The standard deviation for this element is 2014.76 mg ∙ kg⁻¹. Other variables with high standard deviation are phosphorus, iron, magnesium, potassium and manganese. All other parameters have second derivative (SD) values <10.00; cadmium has the lowest standard deviation value of 0.12 mg ∙ kg⁻¹.

Table 3 shows a summary of reflectance data obtained by ADC and the results of all used data normalisation methods. The mean reflectance values for each spectral band are 0.47 for green, 0.41 for red and 0.81 for near-infrared. After MSC data normalisation, the mean red and green band values were slightly changed, while the NIR value was the same in both. The SNV method changed spectral data completely with mean red and green values changed to negative. Application of min–max normalisation effected in NIR spectra became smaller than the red band value. The same was observed for the absorbance values. The standard deviation varied from 0.32 in SNV green band to 0.02 in the NIR MSC band. The MSC method had lower standard deviation values for each band.

Figure 2 presents the relationship between analysed soil characteristics and ADC spectral data in addition to the calculated spectral indices. Correlation values differ from −1.0 (marked as blue on the graph) to 1.0 (red). Most of the soil parameters have a strong negative correlation with spectral data. Only the percentage of sand and some chemical elements, such as Mn, Fe and P, have positive correlation values. Almost every soil parameter is correlated with some of the spectral bands or indices, with the exception of clay and zinc.

Fig. 2

Both regression models were calculated for 21 variables describing soil parameters and 12 variables corresponding to the average reflectance in three bands of the ADC device and the spectral indices calculated on their basis. Additionally, the model was calculated each time for each variant of standardised spectral variables. The data for the Cubist model were divided randomly into a training set of 80% of all data, and a test set which received the remaining 20% of the data. For the PLS model, cross validation type was used which divides the data into segments. The number of segments was set to 10. Regression models, Cubist and PLS, were downloaded in the R software by dedicated packages Cubist, PLS. For the Cubist model, the chosen parameters were the number of committees set to 1 and the number of rules set to 3. The obtained values of the predicted soil parameters were compared with those obtained by laboratory measurements based on the values of the correlation coefficient, root mean square of errors, relative root mean square of errors, RPD and the ratio of yield to interquartile distance (RPIQ) (Table 4).

Yi Peng et al. (2015) used the Cubist model on 328 soil samples to improve SOC modelling at the regional scale. The reference data were a combination of two satellite images and laboratory Vis-NIR measurements. The obtained results were R² = 0.69, RMSE = 2.8, RPD = 1.6 and RPIQ = 0.8. Ng et al. (2022) conducted research that aimed to estimate the available nutrition of many soils using memory-based learning (MBL) algorithm and Cubist regression model. The validation statistics for the prediction of Mehlich III extractable elements using mid-infrared (MIR) spectroscopy are for Ca at R² = 0.91 and RPIQ = 2.11; for Mg R² = 0.82 and RPIQ = 2.39; for K R² = 0.56 and RPIQ = 2.59; for Mn R² = 0.65 and RPIQ = 2.07; for P R² = 0.5 and RPIQ = 1.47; and for Zn R² = 0.59 and RPIQ = 1.59. The result of all elements, except zinc, is similar to or lower than the values obtained in our study. Other studies focused on different methods for soil features estimation. For example, Gholizadeh et al. (2015) focused on using support vector machine regression to establish a relationship between reflectance spectra in the visible near-infrared region and concentrations of Mn, Cu, Cd, Zn and Pb in soil. They used the first and second derivatives (FD and SD), SNV, MSC and continuum removal (CR) to normalise the data. The accuracy was defined by R² and RMSE. FD turned out to be the best normalisation method giving the highest R². For Cu it was 0.78, for Mn 0.6, for Cd 0.8, for Pb 0.68 and for Zn 0.77.

The next step was to answer the question of which soil parameters can be estimated based on multispectral data obtained with ADC and with what accuracy? For that purpose, the threshold values were established for R² ≥0.5 and rel-RMSE ≤0.31. These criteria were met by 12 results (marked in bold in Table 4), which were the contents of sand, silt, SOC, nitrogen, potassium, copper, zinc, lead, manganese and iron, soil pH (pH_H2O, pH_KCL) and CEC. All were obtained using Cubist regression. In Table 4, apart from the regression model, it was also specified which data normalisation method for each parameter guaranteed the highest possible R² value and the lowest possible rel-RMSE value. Normalisation by converting data to absorbance (ABS) most often gave the best results. The MSC method was selected only twice.

According to the model evaluation criteria (Saeys et al. 2005), models for sand, clay, C/N, Pb, Zn and P were considered as not suitable for prediction. The distinction between high and low values was guaranteed by models for K and F. Models for the percentage content of silt, CEC, Cu, Mn, Cd and the first variant of calculating the calcium carbonate content allow for approximate quantitative predictions. According to the given criteria, we can consider pH_H2O, pH_KCL, Mg, Ca and the second method of determining the percentage of calcium carbonate as a good model. Finally, SOC and N were considered as perfect models.

The results, shown in the form of graphs (Fig. 3), present the ratio of the obtained values to the predicted values. The x-axis shows the observed values, i.e. those measured in the laboratory, and the y-axis values are estimated by the regression model.

Fig. 3

For the 12 best-estimated parameters, additional graphs (Fig. 4) were created to illustrate which variables were considered as the most important for building the regression model for each of them. The y-axis contains a list of variables in the order from the one that had the greatest importance for building a given model. The x-axis shows the range of values from 0% to 100%, describing the weight of the variable. Cubist regression model is designed to select only certain parameters for prediction. For this reason, having access to information about VIP can help interpret and understand the given results (Fig. 4).

Fig. 4