Field Proximal Soil Sensor Fusion for Improving High-Resolution Soil Property Maps

Abstract

Mapping soil properties, using geostatistical methods in support of precision agriculture and related activities, requires a large number of samples. To reduce soil sampling and measurement time and cost, a combination of field proximal soil sensors was used to predict and map laboratory-measured soil properties in a 3.4-ha pasture field in southeastern Brazil. Sensor soil properties were measured in situ on a 10 × 10-m dense grid (377 samples) using apparent electrical conductivity meters, apparent magnetic susceptibility meter, gamma-ray spectrometer, water content reflectometer, cone penetrometer, and portable X-ray fluorescence spectrometer (pXRF). Soil samples were collected on a 20 × 20-m thin grid (105 samples) and analyzed in the laboratory for organic C, sum of bases, cation exchange capacity, clay content, soil volumetric moisture, and bulk density. Another 25 samples collected throughout the area were also analyzed for the same soil properties and used for independent validation of models and maps. To test whether the combination of sensors enhances soil property predictions, stepwise multiple linear regression (MLR) models of the laboratory soil properties were derived using individual sensor covariate data versus combined sensor data—except for the pXRF data, which were evaluated separately. Then, to test whether a denser grid sample boosted by sensor-based soil property predictions enhances soil property maps, ordinary kriging of the laboratory-measured soil properties from the thin grid was compared to ordinary kriging of the sensor-based predictions from the dense grid, and ordinary cokriging of the laboratory properties aided by sensor covariate data. The combination of multiple soil sensors improved the MLR predictions for all soil properties relative to single sensors. The pXRF data produced the best MLR predictions for organic C content, clay content, and bulk density, standing out as the best single sensor for soil property prediction, whereas the other sensors combined outperformed the pXRF sensor for the sum of bases, cation exchange capacity, and soil volumetric moisture, based on independent validation. Ordinary kriging of sensor-based predictions outperformed the other interpolation approaches for all soil properties, except organic C content, based on validation results. Thus, combining soil sensors, and using sensor-based soil property predictions to increase the sample size and spatial coverage, leads to more detailed and accurate soil property maps.

1. Introduction

Proximal soil sensing is increasingly used for mapping soils with high spatial detail to support land management and precision agriculture [1,2]. It entails the use of geophysical field and laboratory sensors to measure soil properties (mechanical, electromagnetic, optical, etc.) to directly or indirectly predict and map soil properties of interest. Proximal soil sensors include apparent electrical conductivity (ECa) and magnetic susceptibility (MSa) meters, gamma-ray, X-ray fluorescence and near-infrared spectrometers, and mechanical resistance meters, among others.

Strong correlations have been reported among soil sensor data and soil properties of interest [3,4,5,6]. Usually, these strong correlations are specific to certain soil property-soil sensor combinations. For example, soil salinity is known to increase (or correlate to) soil ECa, thus soil sensors that measure the latter have been used to predict the former [7,8,9]. Going a step further, combining data from different sensors has been shown to assist or improve soil property predictions [4,10,11,12,13,14]. However, the choice of soil sensors to use depends on many factors, including the target soil property, whether they will be used on the field or in the laboratory, the soil and landscape characteristics, cost of the sensor and available budget, and proficiency. Other sensor characteristics like portability, ease of use, measurement support, and the capacity to predict multiple soil properties are also considered, as well as interfering factors and other limitations.

Among the studies that combined proximal sensors for soil property prediction, two in situ ECa sensors (Geonics EM34 and EM38) were used for soil clay prediction in New South Wales, Australia [10]. Gamma-ray variables and ECa, both measured in situ, were combined with elevation, and aerial photos to predict topsoil clay in a 22-ha agricultural field in southwest Sweden, finding a clear superiority of the gamma-ray variables compared to the others [4]. In tropical soils, in situ data from a portable X-ray fluorescence spectrometer (pXRF), a soil color app on a mobile phone, an ECa sensor, and a portable two-band (red, and near-infrared) reflectance sensor were combined to predict many soil chemical and physical properties (K, Mg, Ca, Al, N, C, CEC, particle size fractions, and others) in central Kenya, finding promise in the combination of pXRF and portable reflectance sensors for in situ soil assessments [11]. Furthermore, in southeast Brazil, terrain and parent material variables were combined with magnetic susceptibility and pXRF data, both measured in the laboratory, to predict sand and clay contents in a 150-ha area, and the proximal sensor variables were included in the best models for both soil properties [15].

Albeit previous research has shown the potential of proximal soil sensor combination (fusion) to predict soil properties, a recent review showed that most proximal sensor fusion studies have been done in temperate soils, and that most studies have combined only two or three sensors for soil property prediction and other aims [2]. Thus, the relations among proximally-sensed and laboratory-measured soil properties are still little known when multiple sensors are used in combination, especially when they are used in tropical soils, where only a few proximal soil sensor fusion studies have been done [11,13,15,16], two of which on Amazonian Dark Earths [13,16]. Thus, this study focus on tropical soils and in situ proximal soil sensor fusion aiming to identify which proximal sensors contribute to predict and map chemical and physical soil properties, and how the sensors can be combined to improve the predictions and maps of these properties. The motivation of the study lies in the potential of proximal soil sensors to complement or replace soil sampling and laboratory analyses of these properties, which are costly, time- and energy-consuming, and possibly polluting.

For soil property mapping, ordinary kriging (OK) and its extensions have been widely used in both temperate and tropical regions [17,18,19,20,21,22,23,24,25,26,27,28]. In southeast Brazil, where this study is conducted, OK was used to map ECa measured by a Veris 3100 sensor, as well as corn yield, cost, and profit in a 19-ha irrigated farm [29]. In another study in a 5.7 ha area, OK was compared to inverse distance weighting (IDW) for mapping soil penetration resistance (PR), bulk density (BD), and moisture, using two sampling grids [30]. The best interpolation method varied by soil property and sampling grid. Soil moisture was kriged across a 3.42-ha no-till field (sorghum and soybean) from 102 samples on a 10 × 20 m grid in Reference [31]. However, no cokriging (COK) studies were found in Brazil, which also lacks in situ proximal soil sensing studies relative to other countries [32]. In Brazil, laboratory-based visible and near-infrared reflectance spectroscopy is the most used proximal soil sensing technology, according to Reference [32].

Under this framework, it is hypothesized that: (a) The combination of different proximal soil sensors outperforms individual sensors for soil property prediction. It is assumed that different sensors complement each other in the prediction models, and otherwise, that redundant or useless sensors are not selected in the models; and (b) sensors efficiently improve soil property mapping by indirectly increasing the sample size providing better spatial coverage via soil property prediction. The improved efficiency is achieved by taking sensor measurements directly in situ on a denser sampling grid without collecting, transporting, preparing, and analyzing soil samples in the laboratory.

From the above, the objectives of the study were to:

• Predict six laboratory-measured soil properties from in situ proximal soil sensor covariate data;
• Compare the quality of predictions obtained from individual versus combined proximal soil sensor data;
• Map the six laboratory-measured soil properties using different interpolation approaches;
• Compare the quality of maps derived directly from the observations (raw data) versus those derived from a denser grid of sensor-based predictions.

The six soil properties included three chemical properties—organic C content, sum of bases, and cation exchange capacity—and three physical ones—clay content, volumetric moisture, and bulk density. They encompass a range of soil properties that are important for soil-landscape characterization, understanding soil formation, and guiding land use and management in tropical agriculture.

2. Material and Methods

2.1. Study Area and Soil Sampling

The study was conducted in a 3.4-ha area located in Seropédica, Rio de Janeiro state, southeast Brazil, with central latitude −22.7571 and longitude −43.6630 (Figure 1). The area has a rectangular shape of about 300 m oriented along a toposequence of soils representative of the region, by 140 m across. The toposequence includes Acrisols and Lixisols on the summit and shoulder in the southwest of the area, transitional sandier soils (Acrisols and Planosols) on the slightly undulating backslope in the central portion of the area, and Planosols on the footslope and fluvial terrace in the northeast. The area is under tropical climate, with mean annual temperature and precipitation of 23.2 °C and 1274 mm, respectively [33], and elevations from 22 to 36 m (Figure 1). The land use has been unimproved pasture (Panicum maximum Jacq.) for more than a decade. Soils in the region are formed from granites, gneisses, and migmatites of pre-Cambrian age from the Litoral Fluminense Complex and Serra dos Órgãos Group, with intrusions of basaltic and alkaline rocks of Cretaceous/Tertiary origin, and Quaternary sedimentary deposits from the Piranema Formation [34,35].

A sampling grid of 10 × 10 m (dense grid) was set inside the study area, with 29 points distributed along the toposequence by 13 points across (377 points), leaving a 10-m buffer outside the grid (Figure 1). Six proximal soil sensors (Figure 2) were used to take in situ measurements on these 377 sites. Apparent magnetic susceptibility (MSa) [36] was measured by the KT-10 S/C sensor (Terraplus Inc., Richmond Hill, ON, Canada; Figure 2A, sensor #4). Apparent electrical conductivity (ECa) [37] was measured at the surface by the KT-10 S/C sensor, and at 0–10, 0–20, and 0–40 cm by the Rabellis sensor (Embrapa Instrumentação Agropecuária, São Carlos, SP, Brazil; Figure 2B) [38]. The RS-230 BGO gamma-ray spectrometer (Radiation Solutions Inc., Mississauga, ON, Canada; Figure 2A, sensor #5) [39] was used to measure the dose rate and equivalent of U (eU) and Th (eTh) contents, with measurements taken over 120 s at the soil surface. The soil volumetric moisture (θ) was measured by the CS650 water content reflectometer (WCR) (Campbell Scientific Inc., Logan, UT, USA; Figure 2A, sensor #3) [40], and the cone penetration resistance (PR) by the PenetroLOG cone penetrometer (Falker Automação Agrícola Ltd.a., Porto Alegre, RS, Brazil; Figure 2A, sensor #1), using a type 2 cone, both at 0–10 cm. Finally, the DP-6000 pXRF sensor (Olympus Scientific Solutions Americas Inc., Waltham, MA, USA; Figure 2A, sensor #2) [41] was used in “Soil” mode to measure total element contents of various elements at the soil surface. The three beams were used for 30 s each, totaling 90 s per sample. Many element contents were below the sensor detection limits for many samples. Thus, only those element contents with at least 374 valid readings were selected for the study, including K, Ti, Mn, Fe, Zn, Rb, Sr, Zr, Ba, Cr, and Pb. In previous studies, this pXRF sensor achieved 72–90% recovery against NIST certified soils for Fe, 85–95% for Mn [42], and 90–109% for Ti [43].

From the dense grid, a grid of 20 × 20 m (thin grid), with 15 by 7 points (105 points), was set by skipping every other sampling row and column (Figure 1). Additionally, 25 sampling sites for independent validation were distributed across the study area by conditioned Latin hypercube sampling (cLHS) [44], using the sampling row, sampling column, and elevation as strata. On these 130 sites, soil samples were collected at 0–10 cm and analyzed in the laboratory, according to Reference [45], for organic C (OC) content by wet oxidation with K₂Cr₂O₇ + H₂SO₄, exchangeable bases (Ca and Mg extracted with KCl, and K and Na extracted with HCl + H₂SO₄), exchangeable acidity (H and Al extracted with Ca(C₂H₃O₂)₂), clay content by the hydrometer method, and soil θ and BD measured from 100-cm³ steel-ring samples. The sum of bases (SB) was derived as the sum of exchangeable bases, and the cation exchange capacity (CEC) as the sum of SB and exchangeable acidity.

2.2. Predictive Modeling and Mapping

The soil chemical (OC, SB, and CEC) and physical properties (clay, θ, and BD) were modeled by multiple linear regression (MLR) as a function of data from one specific sensor, and combined data from multiple sensors, respectively. Data were log-transformed when positively skewed, and stepwise variable selection with p < 0.05 was applied. Descriptive statistics of individual variables, and linear correlation coefficients (r in Equation (1)) among individual variables were calculated. The 105 samples on the thin grid were used to train the models, whereas the 25 cLHS samples were used as independent validation samples to test and compare results on an external data set. The sensor-measured θ has many missing values from the first field campaign; thus, it was included separately in the models for comparison only, and not used to produce soil property maps.

Exceptionally, the pXRF sensor was evaluated alone and not combined with the other sensors. The pXRF sensor alone provided many predictor variables, as it measures the contents of multiple elements simultaneously, and some of these element contents showed moderate to strong correlations with some target soil properties. As such, it was anticipated that the pXRF alone could outperform the other sensors individually, and even the other sensors combined, for predicting the soil properties, making it reasonable to evaluate it separately. The adjusted coefficient of determination (R²_adj in Equation (2)) was used to assess model fit, and the root mean square error (RMSE in Equation (3)) to assess prediction accuracy. The smallest RMSE calculated on the 25 cLHS validation samples was used to select the best model for each laboratory soil property, respectively.

$$r={\displaystyle \sum }_{i=1}^N\left(x_i-\overline{x}\right)\left(y_i-\overline{y}\right)/\sqrt{{\displaystyle \sum }_{i=1}^N{\left(x_i-\overline{x}\right)}^2{\displaystyle \sum }_{i=1}^N{\left(y_i-\overline{y}\right)}^2}$$

(1)

$$R^2{}_{adj}=1-\left[\left(1-r^2\right)\left(N-1\right)/\left(N-p-1\right)\right]$$

(2)

$$RMSE=\sqrt{{\displaystyle \sum }_{i=1}^N{\left(y_i-{\widehat{y}}_i\right)}^2/N}$$

(3)

where $r$ is the linear correlation coefficient, $R^2{}_{adj}$ is the adjusted coefficient of determination, $RMSE$ is the root mean square error, $x_i$ and $y_i$ are the observed values of $x$ and $y$ variables, respectively, $\overline{x}$ and $\overline{y}$ are their respective mean values, $N$ is the sample size, $p$ is the number of predictors in the model, and ${\widehat{y}}_i$ are the predicted values of $y$.

Three interpolation approaches were compared to map the six above mentioned soil properties with 1-m spatial resolution. The first approach was to interpolate the 105 observations (raw data) from the thin grid using OK [46]. These were considered baseline maps—that is, maps that are produced without the use of ancillary proximal sensor data. In the second approach, first predictions were made for the six laboratory-measured soil properties from their best prediction models, respectively, on the 352 dense grid sites (377 sites minus the 25 cLHS sites set apart for validation; empty circles in Figure 1), which contain sensor data only. Then, the 352 predictions on the dense grid were interpolated by OK. From the second approach, two maps were derived for each soil property, one from the best single-sensor model, and another from the combined-sensor model. The third approach was to use ordinary COK [46] to interpolate the observations from the thin grid aided by sensor covariate data from the dense grid. For the latter, the proximally-sensed property that had the highest correlation and an adequate spatial cross-correlation structure with the target soil property was selected as the covariate. The empirical variograms (Equation (4) [46]) of all soil properties were fit by the spherical variogram model (Equation (5) [46]) using ordinary least squares or manually. The quality of the predictions from the three interpolation approaches was assessed by calculating the RMSE on the 25 validation samples and then compared.

$$\gamma =\left(1/2N_h\right){\displaystyle \sum }_{i=1}^{N_h}{\left[z\left({\alpha }_i\right)-z\left({\alpha }_i+h\right)\right]}^2$$

(4)

$$\widehat{\gamma }=\{\begin{array}{c}{c}_0+c\left[1.5\left(h/a\right)-0.5{\left(h/a\right)}^3\right]\\ c_0+c\end{array}\begin{array}{c}, if h\le a\\ , if h>a\end{array}$$

(5)

where $\gamma$ is the observed semivariance of variable $z$ at lag distance $h$, $z\left({\alpha }_i\right)$ and $z\left({\alpha }_i+h\right)$ are the observed values of $z$ at $N_h$ pairs of locations separated by lag distance $h$, and $\widehat{\gamma }$ is the fitted spherical variogram of $z$ as a function of lag distance $h$ and variogram parameters nugget effect ($c_0$), sill ($c$), and range ($a$).

3. Results and Discussion

3.1. Descriptive Statistics and Correlations

Among the chemical soil properties, OC varied from 3.7 to 28 g kg⁻¹, with a mean of 11.3 g kg⁻¹, whereas the CEC varied from 2.6 to 12.9 cmol_c kg⁻¹, with a mean of 6.8 cmol_c kg⁻¹ (

Table 1. Descriptive statistics of laboratory-measured, and field proximally-sensed soil properties.

Property 1	N 1	Min 1	Max 1	Mean	Median	SD 1	Skew 1
OC (g kg−1)	130	3.7	28.0	11.3	11.4	3.9	0.64
SB (cmolc kg−1)	130	1.4	8.4	3.9	3.8	1.5	0.46
CEC (cmolc kg−1)	130	2.6	12.9	6.8	6.8	1.9	0.29
Clay (g kg−1)	130	20	380	176	160	93	0.24
θ (%)	130	4.2	31.7	15.4	15.6	5.6	0.25
BD (g cm−3)	130	1.12	1.72	1.52	1.54	0.09	−1.18
WCR θ (%)	310	2.3	35.7	14.2	14.1	6.7	0.42
KT MSa (10−3 SI)	376	0.0	3.8	0.3	0.2	0.3	4.44
KT log(MSa)	376	−3.9	1.3	−1.4	−1.5	0.8	−0.06
KT ECa (S m−1)	376	0.0	6.2	1.5	1.4	1.0	0.94
Rab ECa 0–10 (S m−1)	374	0.0	52.7	1.9	0.3	4.5	5.99
Rab log(ECa 0–10)	374	−5.1	4.0	−0.8	−1.3	1.6	0.55
Rab ECa 0–20 (S m−1)	374	0.1	99.9	6.7	3.8	11.5	4.50
Rab log(ECa 0–20)	374	−2.9	4.6	1.2	1.3	1.2	−0.36
Rab ECa 0–40 (S m−1)	374	0.0	10.7	1.4	1.1	1.4	3.26
Rab log(ECa 0–40)	374	−3.0	2.4	0.0	0.1	0.9	−0.30
RS DR (μR h−1)	376	0.2	5.2	2.7	2.7	0.9	0.33
RS eU (mg kg−1)	376	0.3	3.4	1.3	1.2	0.4	0.59
RS eTh (mg kg−1)	376	0.0	13.9	6.2	5.9	2.6	0.36
PR (kPa)	376	53	2697	810	678	493	1.28
pXRF K (mg kg−1)	376	0	5302	1997	1866	883	0.81
pXRF Ti (mg kg−1)	376	2254	9638	5141	4947	1465	0.66
pXRF Mn (mg kg−1)	376	83	3421	506	401	375	3.17
pXRF Fe (mg kg−1)	376	2794	53,946	14,014	12,251	8193	1.48
pXRF Zn (mg kg−1)	376	5	153	40	35	24	1.73
pXRF Rb (mg kg−1)	376	3	65	22	21	12	0.89
pXRF Sr (mg kg−1)	376	5	66	25	23	11	0.73
pXRF Zr (mg kg−1)	376	319	2251	917	857	338	1.17
pXRF Ba (mg kg−1)	376	69	503	216	203	78	0.86
pXRF Cr (mg kg−1)	375	0	40	10	9	8	0.94
pXRF Pb (mg kg−1)	374	3	34	14	14	5	0.56

¹ N, number of observations; Min, minimum; Max, maximum; SD, standard deviation; Skew, skewness coefficient; OC, organic C content; SB, sum of bases; CEC, cation exchange capacity; θ, volumetric moisture; BD, bulk density; MSa, apparent magnetic susceptibility; ECa, apparent electrical conductivity; WCR, CS650 water content reflectometer; KT, KT-10 S/C MSa/ECa sensor; Rab, Rabellis ECa sensor; 0–10…40, depths of measurement; RS, RS-230 BGO gamma-ray spectrometer; DR, dose rate; eU, equivalent U content; eTh, equivalent Th content; PR, mean cone penetration resistance at 0–10 cm; pXRF, DP-6000 portable X-ray fluorescence spectrometer.

). Among the physical properties, the clay content varied between 20 and 380 g kg⁻¹, and soil θ between 4.2 and 31.7%. These values agree to those reported previously for Planosols [47] and Acrisols [48] of the region.

Most laboratory- and sensor-measured soil properties were significantly correlated (

Table 2. Linear correlations among laboratory- and sensor-measured soil properties.

Property 1	OC	SB	CEC	Clay	θ	BD
OC	1
SB	0.67 *	1
CEC	0.76 *	0.87 *	1
Clay	0.65 *	0.68 *	0.62 *	1
θ	0.63 *	0.63 *	0.62 *	0.85 *	1
BD	−0.45 *	−0.22 *	−0.31 *	−0.20 *	−0.13ns	1
WCR θ	0.63 *	0.56 *	0.54 *	0.68 *	0.76 *	−0.12ns
KT log(MSa)	0.34 *	0.50 *	0.42 *	0.53 *	0.50 *	−0.09ns
KT ECa	0.38 *	0.41 *	0.43 *	0.45 *	0.49 *	−0.16ns
Rab log(ECa 0–10)	0.45 *	0.47 *	0.49 *	0.35 *	0.47 *	−0.01ns
Rab log(ECa 0–20)	0.49 *	0.51 *	0.54 *	0.38 *	0.49 *	−0.05ns
Rab log(ECa 0–40)	0.49 *	0.54 *	0.54 *	0.43 *	0.59 *	−0.02ns
RS DR	0.51 *	0.51 *	0.52 *	0.70 *	0.60 *	−0.13ns
RS eU	−0.01ns	−0.10ns	0.02ns	−0.09ns	−0.09ns	−0.02ns
RS eTh	0.54 *	0.57 *	0.54 *	0.78 *	0.67 *	−0.12ns
PR	0.32 *	0.36 *	0.28 *	0.50 *	0.44 *	0.08ns
pXRF K	0.28 *	0.33 *	0.33 *	0.34 *	0.22 *	−0.15ns
pXRF Ti	−0.11ns	−0.16ns	−0.08ns	−0.13ns	−0.16ns	0.06ns
pXRF Mn	0.13ns	0.21 *	0.27 *	0.00ns	0.00ns	−0.05ns
pXRF Fe	0.43 *	0.49 *	0.38 *	0.80 *	0.67 *	−0.08ns
pXRF Zn	0.24 *	0.38 *	0.29 *	0.33 *	0.23 *	−0.07ns
pXRF Rb	0.44 *	0.56 *	0.49 *	0.69 *	0.54 *	−0.17 *
pXRF Sr	0.47 *	0.46 *	0.49 *	0.29 *	0.26 *	−0.24 *
pXRF Zr	−0.33 *	−0.49 *	−0.36 *	−0.53 *	−0.50 *	0.03ns
pXRF Ba	0.46 *	0.46 *	0.43 *	0.72 *	0.57 *	−0.14ns
pXRF Cr	0.20 *	0.29 *	0.27 *	0.37*	0.22 *	−0.19 *
pXRF Pb	0.62 *	0.43 *	0.50 *	0.60 *	0.46 *	−0.30 *

¹ OC, organic C content; SB, sum of bases; CEC, cation exchange capacity; θ, volumetric moisture; BD, bulk density; MSa, apparent magnetic susceptibility; ECa, apparent electrical conductivity; WCR, CS650 water content reflectometer; KT, KT-10 S/C MSa/ECa sensor; Rab, Rabellis ECa sensor; 0–10…40, depths of measurement; RS, RS-230 BGO gamma-ray spectrometer; DR, dose rate; eU, equivalent U content; eTh, equivalent Th content; PR, mean cone penetration resistance at 0–10 cm; pXRF, DP-6000 portable X-ray fluorescence spectrometer. *, significant at the 0.05 significance level; ns, not significant.

). The highest correlations among laboratory-measured properties were found between SB and CEC (r = 0.87), clay content and soil θ (0.85), and OC and CEC (0.76). The sensor variables with the highest correlations with laboratory-measured properties were eTh from the gamma-ray sensor with clay content (0.78), and θ from the water content sensor with laboratory-measured θ (0.76). The sensor-measured properties θ, dose rate and eTh had moderate to high correlations (r > 0.50) with all laboratory-measured properties except for BD, indicating potential variables for proximal sensor fusion. This potential was confirmed by their inclusion in the combined-sensor models, as discussed in the next section. In detail, the sensor-measured θ had the highest correlation among sensor variables with laboratory-measured OC and θ, and second-highest with SB (after eTh); eTh had the highest correlation with SB, and second-highest with clay (after pXRF Fe) and θ (after WCR θ); and both had the highest correlation with CEC (Table 2).

3.2. Individual-Versus Combined-Sensor Models

The MLR prediction models had moderate to good fits with R²_adj > 0.50 for all soil properties except for BD (

Table 3. Model training and validation results for the laboratory soil properties predicted as a function of data from individual and combined proximal soil sensors.

Property 1	Model 1	Selected Covariates 1	Training 1			Validation 1
			Nt	R2adj	RMSEt	Nv	RMSEv 2
Individual sensors (excluding WCR)
OC (g kg−1)	pXRF	K, Ti, Mn, Zn, Rb, Sr, Zr, Ba, Cr, Pb	102	0.58	2.2	25	2.5
	KT	log(MSa), ECa	104	0.11	3.3	25	3.1
	RS	DR, eTh	103	0.33	2.8	25	3.2
	Rab	log(ECa 0–10), log(ECa 0–40)	103	0.30	2.8	24	3.3
SB (cmolc kg−1)	pXRF	Ti, Fe, Rb, Sr, Zr	99	0.64	0.8	25	0.9
	RS	DR, eTh	105	0.27	1.2	25	0.9
	Rab	log(ECa 0–10), log(ECa 0–40)	105	0.25	1.3	24	1.1
	KT	log(MSa), ECa	103	0.24	1.6	25	2.9
CEC (cmolc kg−1)	RS	DR, eTh	105	0.27	1.7	25	1.4
	KT	log(MSa), ECa	103	0.24	1.6	25	1.5
	Rab	log(ECa 0–10), log(ECa 0–40)	105	0.27	1.7	24	1.5
	pXRF	Ti, Sr, Zr, Pb	102	0.52	1.3	25	1.7
Clay (g kg−1)	pXRF	Ti, Fe, Zn, Rb, Zr, Cr, Pb	103	0.81	37	25	40
	RS	eU, eTh	100	0.79	38	25	60
	KT	log(MSa), ECa	105	0.21	78	25	80
	Rab	log(ECa 0–10), log(ECa 0–40)	105	0.16	81	24	97
θ (% m/v)	pXRF	Ti, Zn, Zr, Ba, Cr, Pb	103	0.57	3.2	25	5.2
	RS	DR, eTh	105	0.44	3.8	25	5.4
	Rab	log(ECa 0–10), log(ECa 0–40)	105	0.34	4.1	24	6.0
	KT	log(MSa), ECa	105	0.22	4.4	25	6.1
BD (g cm−3)	pXRF	K, Ti, Mn, Fe, Zn, Rb, Sr, Zr, Cr, Pb	102	0.01	0.078	25	0.085
	RS	DR, eTh	104	0.00	0.081	25	0.090
	Rab	log(ECa 0–10), log(ECa 0–40)	103	−0.02	0.077	24	0.091
	KT	log(MSa), ECa	104	0.02	0.081	25	0.093
Combined sensors (excluding pXRF)
OC (g kg−1)	Combined	log(Rab ECa 0–20), eTh	103	0.46	2.5	24	3.0
	Combined + WCR	WCR θ, DR, PR	46	0.58	2.2	24	3.3
SB (cmolc kg−1)	Combined	log(MSa), log(Rab ECa 0–40), DR	105	0.39	1.1	24	0.7
	Combined + WCR	WCR θ, log(MSa)	47	0.28	1.2	24	0.7
CEC (cmolc kg−1)	Combined	ECa, log(Rab ECa 0–20), DR	104	0.47	1.4	24	1.2
	Combined + WCR	WCR θ, DR	46	0.42	1.4	24	1.6
Clay (g kg−1)	Combined + WCR	WCR θ, log(MSa), log(Rab ECa 0–20), DR, eU, PR	44	0.94	21	23	47
	Combined	KT ECa, eU, eTh, PR	101	0.79	38	25	54
θ (% m/v)	Combined + WCR	WCR θ, log(ECa 0–20), log(ECa 0–40), DR	47	0.81	2.2	23	3.8
	Combined	log(Rab ECa 0–20), log(Rab ECa 0–40), eTh	104	0.64	3.0	24	4.7
BD (g cm−3)	Combined	eTh, PR	103	0.07	0.074	25	0.090
	Combined + WCR	WCR θ, eTh	46	0.10	0.066	24	0.095

¹ Nt, number of observations in the training set; R²_adj, adjusted coefficient of determination; RMSEt, root mean square error of training; Nv, number of observations in the validation set; RMSEv, root mean square error of validation; OC, organic C content; SB, sum of bases; CEC, cation exchange capacity; θ, volumetric moisture; BD, bulk density; MSa, apparent magnetic susceptibility; ECa, apparent electrical conductivity; KT, KT-10 S/C MSa/ECa sensor; Rab, Rabellis ECa sensor; RS, RS-230 BGO gamma-ray spectrometer; pXRF, DP-6000 portable X-ray fluorescence spectrometer; WCR, CS650 water content reflectometer; 0–10…40, depths of measurement; DR, dose rate; eTh, equivalent Th content; eU, equivalent U content; PR, mean cone penetration resistance at 0–10 cm. ² For each target soil property, the RMSEv is ordered from the smallest (best) to the largest (worst) in the “Individual sensors” and “Combined sensors” sections, respectively.

). The highest R²_adj were found for clay content, and soil θ, both as a function of combined sensors plus the CS650 WCR water content sensor, with R²_adj of 0.94 and 0.81, respectively, stressing the importance of including a water content sensor in proximal sensor combinations, as previously noted by References [3,49]. However, clay content and θ were the only soil properties that benefited from adding WCR θ to the combined models. Considering the RMSE of external validation to select the best models, the pXRF covariates (element contents) derived the best models for OC, clay, and BD, whereas combined models (with or without WCR θ) outperformed the individual sensors, including pXRF, for SB, CEC, and θ. Among individual sensors, the pXRF, gamma-ray, and water content sensors showed the best performances for soil property assessment.

The pXRF sensor was superior to the other sensors combined in modeling three out of six soil properties, and was the best among individual sensors to predict all soil properties except for CEC, based on the RMSE of validation (shown in Table 3 in order from the best to the worst for each target property, respectively). A possible explanation is the fact that this sensor measures the contents of many chemical elements simultaneously, and thus, provides many covariates for soil property prediction. Moreover, some correlations among pXRF covariates and target soil properties were moderate to strong, for example, between OC and Pb (r = 0.62); SB and Rb (0.56); clay and Fe (0.80), Ba (0.72), Rb (0.69) and Pb (0.60); and θ and Fe (0.67), Ba (0.57) and Rb (0.54) (Table 2). On the other hand, the relatively lower correlations observed among the target soil properties and pXRF K, Ti and Mn were not expected, since these elements constitute more common soil minerals and are more abundant in soils than Ba, Rb and Pb, for example (Table 1). The gamma-ray sensor was the second-best among individual sensors for predicting all soil properties, and the best one for CEC prediction (Table 3), which is supported by the relatively high correlations observed between the sensor-measured eTh and dose rate, and all target soil properties except for BD (Table 2).

All pXRF models selected Ti and Zr as covariates (Table 3). Titanium is the second-most abundant element among those measured by pXRF (Table 1), followed by K, which was included in the OC and BD models only, and then Zr. Titanium and Zr had a moderately high correlation (r = 0.63), whereas K had its highest correlations with Rb (0.78) and Sr (0.64). Titanium and Zr are usually present in primary minerals (rutile—TiO₂, ilmenite—FeTiO₃, and zircon—ZrSiO₄) in the sand-size fraction of soils [50]. Thus, they can indicate different parent materials, lithological discontinuities, and/or preferential weathering responding to water and relief dynamics. In the study area, Ti and Zr had 0.33 and 0.60 correlations with fine sand content (0.05–0.20 mm), respectively, with both Ti and Zr, and fine sand content presenting higher values in the fluvial terrace and lower values on the summit. Presumably, the sand on the terrace present at the topsoil (A + E horizons) derive from other parent materials carried in from outside the study area by turbulent water flows, which corroborates the correlated Ti and Zr patterns. Unexpectedly, Fe—the most abundant among the measured elements and notably important in tropical soils forming Fe oxides—was not selected in the OC, CEC, and soil θ MLR models, although its correlation to these properties was higher compared to some of the model-selected pXRF elements. Iron was strongly correlated with Ba (0.89), which may explain why they were mutually exclusive in the models. It should be mentioned that the three pXRF models that outperformed the combined models (for OC, clay, and BD) were the ones that selected the most pXRF elements, up to 10 out of 11 elements (Table 3).

Portable X-ray fluorescence has been used to predict sand, and clay contents in Louisiana, USA, obtaining R² of 0.86 and 0.96, respectively [51]. Moreover, in Louisiana, USA, soil pH was predicted from pXRF data, obtaining an R² of 0.77 [5], and in California and Nebraska, USA, CEC was predicted with R² of 0.91 [6]. In the aforementioned studies, samples were measured in the laboratory using a pXRF sensor, whereas in the present study pXRF measurements were taken in situ. Comparatively, the R²_adj of the pXRF models for clay and CEC obtained in this study were 0.81 and 0.52, respectively. When used to predict the same pXRF elements, but measured by inductively coupled plasma optical emission spectrometry (ICP-OES), the field-measured pXRF elements predicted Fe, Mn, Zn, and Pb contents with R² of up to 0.91, 0.70, 0.61 and 0.66, respectively [52]. Note that the pXRF sensor measures many other soil elements that were not included in this study, either because they were below the detection limit or because they are not provided in “Soil” mode, which was used.

All combined models without the CS650 WCR water content sensor selected at least one property measured by the gamma-ray sensor, with a preference for dose rate or eTh, with eU only selected for clay (Table 3). All combined + WCR models selected WCR θ in the stepwise process. Apparent electrical conductivity variables were included in all combined models without WCR; however, when WCR θ was included, the ECa variables were left out, except for clay and θ. Penetration resistance was included in the prediction models for OC, clay, and BD. The combination of proximal soil sensors in these models improved predictions (that is, reduced the RMSE of validation) of all soil properties investigated compared to individual sensors, disregarding the pXRF sensor (which was not considered for this comparison because it measures many soil elements at once). These results show that proximal sensor fusion is superior to individual sensors for soil property prediction. As an exception, BD could not be well predicted from any sensor or sensor combination, as its validation RMSE was as high as its standard deviation (BD varies little in the study area) (Table 1 and Table 3). In other words, its ratio of validation RMSE to standard deviation, commonly used to infer prediction quality, was close to 1, indicating weak prediction capacity. For all other soil properties, the validation RMSE was lower than their standard deviation, respectively.

The combination of sensors improving soil property predictions relative to single sensors was also observed in Reference [4], who obtained better clay content predictions by combining ECa and gamma-ray sensors. Their modeling efficiency (similar to R²) improved from 0.85 and 0.94 (ECa-only and gamma-ray-only models, respectively) to 0.96 when the sensors were combined. The ECa and gamma-ray sensor combination also better-discriminated soils that had similar outputs to one individual sensor [53]. In Brazilian Latosols, R²_adj of 0.67 and 0.58 were obtained for clay and sand content prediction, respectively, by combining magnetic susceptibility and pXRF data [15]. Besides improving soil property predictions, the combination of proximal soil sensors offers other advantages, due to the complementary features of individual sensors. For instance, PR required a water content sensor for its proper interpretation as a measure of soil compaction [54], or for its prediction from a mechanical resistance sensor [3].

3.3. Baseline Versus Sensor-Aided Maps

Compared to OK from the 105 thin grid raw observations (that is, the baseline maps), the interpolated sensor-based predictions on the dense grid with 352 sites achieved the least validation errors, except for OC (

Table 4. Fitted spherical variogram parameters of laboratory soil properties, and external validation errors of produced maps.

Property 1	Method 1	Grid	Model/Covariate 1	Nugget	Sill	Range (m)	Nugget/Sill (%)	RMSEv 1,2
OC (g kg−1)	OK	Thin		5.00	28.00	348	17.9	2.8
	OK	Dense	pXRF	3.16	12.07	291	26.2	2.9
	OK	Dense	Combined	2.50	7.50	269	33.3	3.0
	COK	Mixed	RS eTh	8.33	17.57	220	47.4	3.2
SB (cmolc kg−1)	OK	Dense	Combined	0.30	1.20	225	25.0	0.9
	OK	Thin		1.05	3.15	254	33.3	1.0
	OK	Dense	pXRF	0.31	1.73	147	18.0	1.0
	COK	Mixed	RS eTh	1.70	2.79	230	61.0	1.3
CEC (cmolc kg−1)	OK	Dense	Combined	0.80	2.30	138	34.8	1.5
	OK	Dense	RS	0.32	1.58	281	20.3	1.5
	OK	Thin		2.00	6.50	339	30.8	1.6
	COK	Mixed	RS eTh	2.67	4.50	220	59.2	1.9
Clay (g kg−1)	OK	Dense	pXRF	925	11573	228	8.0	57
	OK	Thin		450	13915	278	3.2	60
	COK	Mixed	RS eTh	2408	12412	220	19.4	61
	OK	Dense	Combined	1118	9267	294	12.1	63
θ (% m/v)	OK	Dense	Combined	7.00	32.00	285	21.9	5.6
	OK	Thin		9.33	38.39	275	24.3	5.8
	COK	Mixed	RS eTh	10.58	31.97	205	33.1	6.1
	OK	Dense	pXRF	4.72	19.22	190	24.6	6.3
BD (g cm−3)	OK	Dense	pXRF	0.0002	0.0007	78	31.7	0.086
	OK	Dense	Combined	0.0004	0.0005	112	74.9	0.087
	OK	Thin		0.0065	0.0100	316	65.0	0.089
	COK	Mixed	KT ECa	0.1971	0.1988	245	99.1	0.210

¹ RMSEv, root mean square error of validation; OC, organic C content; SB, sum of bases; CEC, cation exchange capacity; θ, volumetric moisture; BD, bulk density; eTh, equivalent Th content; ECa, apparent electrical conductivity; OK, ordinary kriging; COK, ordinary cokriging; pXRF, DP-6000 portable X-ray fluorescence spectrometer; RS, RS-230 BGO gamma-ray spectrometer; KT, KT-10 S/C MSa/ECa sensor. ² For each target soil property, the RMSEv is ordered from the smallest (best) to the largest (worst).

), showing that, overall, proximal soil sensors produces better soil property maps. In detail, soil property predictions from the combined models derived the most accurate maps for soil SB, CEC, and θ, whereas clay, and BD maps benefited most from the pXRF predictions, according to the smaller RMSE of validation (Table 4). Only in two cases, COK performed a little better than any other approach: It was better than OK from combined sensor-based predictions for clay content; and better than OK from pXRF-based predictions for soil θ. Equivalent Th, measured by the gamma-ray sensor, was selected as the covariate for COK of all properties except BD (which used KT-10 ECa) because it had a high correlation and an adequate spatial cross-dependence structure with them. These results stress the superiority of OK of sensor-based predictions over the other approaches.

The empirical and fitted variograms of the soil properties are shown in Figure 3. From left to right, the variograms were derived, respectively, from: The 105 thin grid raw samples; the 352 dense grid predictions from the best individual-sensor model; the 352 dense grid predictions from the combined-sensor model; and the 105 thin grid samples and collocated dense grid samples of the selected covariate. The best prediction models, the selected COK covariates, and the variogram parameters of the fitted spherical models are listed in Table 4. For each target soil property, similar spatial dependence structures were found among the three interpolation approaches (Table 4; Figure 3). However, in general, the variogram range decreased when data from the dense grid was considered, either using OK or COK.

Among OK approaches, except for clay content, the nugget variance decreased when the dense grid sensor data was used, which makes sense because data is available to explain the variance over shorter distances. Most soil properties showed a marked increase in the empirical semivariance at around 200 m, which is related to the approximate distance between the hillslope positions with the most distinct soil properties values in the area. In detail, on the summit and shoulder, in the southwest of the area (Figure 1), Lixisols are present with higher OC, SB, CEC, and θ values, whereas on the footslope and floodplain, in the central-east part, Planosols with lower soil property values are found. This increase in the semivariance can indicate an uncorrected geographic trend, which is not desired for OK. Thus, universal kriging using the geographic coordinates and elevation to remove the trend was tested in a separate study, and did not show considerable differences in the spatial patterns and quality of output maps that supported its use instead of OK [55].

The derived soil property maps are shown in Figure 4. From left to right, the maps were produced by: OK of the 105 thin grid samples; OK of the 352 dense grid predictions from the best individual-sensor model; OK of the 352 dense grid predictions from the combined-sensor model; and COK of the 105 thin grid samples aided by the 352 dense grid samples. The best prediction models and the selected COK covariates are listed in Table 4. The maps show similar spatial trends among soil properties, with OC, SB, CEC, clay, and θ showing trends in the same direction, and BD showing opposite ones (Figure 4). One of the reasons is their significant correlations (Table 2) and another one is their similar spatial dependence structures (Table 4; Figure 3). These correlations and the similarities among soil property maps derive from the common soil-forming factors and processes acting in the study area, notably those related to the interconnected water dynamics and relief processes and patterns that resulted in the different geomorphic/landscape strata in the area.

In detail, the upslope portion (summit and shoulder) in the southwest of the area has the highest OC, SB, CEC, clay, and θ values, the latter due to the higher clay content (Figure 4). In the east, there is a floodplain with a drainage channel cutting across it in the south to north direction (Figure 1), with sand-rich deposits (fluvial terrace) on both sides of the channel where the lowest values of these soil properties are found. On the backslope in the central portion, the soil property values are intermediate. Due to the combined effect of turbulent water flows with high energy in the events of heavy rains, seasonal flooding, and fluctuating water table in a depositional environment, the Planosols on the fluvial terrace undergo processes of preferential deposition of sand-sized particles from outside areas, and removal of clays, which are washed out by the drainage water, forming E horizons. The channel drains outward to the extreme north of the area, where high values are found again for these soil properties. Bulk density spatial patterns were the opposite of clay (and thus, of the other properties), as sandier soils have higher BD. Thus, high BD values concentrated on the footslope and sand-rich terrace, whereas low BD values occurred on the summit and extreme north of the area (Figure 4).

For each soil property, the quality of predictions was reasonably similar among interpolation approaches (Table 4), meaning that only modest improvement in map accuracy was achieved by using sensor covariates. However, the sensor-aided derived maps were more detailed (Figure 4), since the sensor-based predictions had better spatial coverage with more samples to be interpolated. On the other hand, COK did not perform as well as OK of sensor-based predictions, as similarly found by References [17,56]. In Reference [17], kriging predictions from terrain-based models using two approaches (regression-kriging models “A” and “B”) were compared to three other methods, including COK, to predict the depth of solum, depth to bedrock, topsoil gravel, and subsoil clay content in a sub-catchment area in Adelaide, Australia. For all properties, kriging the terrain-based predictions using either regression-kriging model outperformed COK and three other methods tested. Sensor-derived ECa data was used to predict soil horizon depth in a 97-ha grassland area in the Netherlands [56]. The predictions obtained by OK were better than those from COK of the horizon depth aided by ECa data. In contrast, the best predictions for sand, silt, and clay contents were obtained by COK using ECa as a covariate, relative to three univariate interpolation methods in 10 to 45-ha farm fields in Poland [23]. However, they did not interpolate sensor-based predictions.

Different sampling designs were also compared for mapping various soil properties in References [30,57]. In Reference [30], similar to our study, soil PR, BD, and soil moisture maps derived from a dense grid of 145 samples (20 × 20 m) were more accurate than their corresponding maps derived from a thin grid of 41 samples (40 × 40 m). When using the sensor-based prediction on the dense grid for OK, the RMSE of validation improved, indicating more accurate maps. In Reference [57], grid sampling was compared to simple random sampling for kriging nine soil properties in a 10-ha tobacco field in southwest China. Using 115 samples from each design, grid sampling provided more accurate maps of soil organic matter, CEC, total N, total K and available K, simple random sampling was preferred for soil pH, total P, and available P, and no difference among designs was found for nitrate N.

3.4. Recommendations and Final Considerations

For individual use, among the sensors tested, the DP-6000 pXRF and the RS-230 BGO gamma-ray sensors are recommended for their portability, ease of use, and potential to predict and aid fine-scale mapping of other soil properties. Moreover, the pXRF sensor derived the best predictive models for OC, clay, and BD, the eTh (measured by the gamma-ray sensor) had relatively high correlations with all soil properties except BD. It is also recommended to include the CS650 WCR or another water content sensor to improve soil property predictions. In this study, the WCR θ had the highest correlations with OC, CEC, and soil θ among proximal sensor variables (Table 2), and was selected as a covariate in all combined + WCR stepwise models (Table 3). A water content sensor is important for interpreting and/or correcting other sensor measures that are influenced by soil moisture (for example, PR and ECa).

The pXRF, gamma-ray, and other sensors used in the study display the measurements as they are taken in the field. As such, beyond soil property prediction and mapping, potential applications for these sensors include: Rapid in situ soil characterization; assessment of the horizontal and vertical variation (or homogeneity) of soils; delineation of soil mapping units or management zones; assessment of soil chemical and physical quality and limitations; allocation of samples for soil mapping; identification of soil and mineral anomalies in geochemical surveys; and others.

Via MLR prediction, the proximal soil sensors indirectly increased the sample size improving the spatial coverage for interpolation, which provided more accurate and more detailed soil property maps. Thus, for soil property mapping, proximal soil sensors could be used in place of laboratory soil analysis, especially when used in combination, in similar soil-landscape settings. This would save time, personal, financial, and environmental resources, as sample collection, transportation, preparation, and laboratory analysis would be avoided. Moreover, it is expected that this successful application of proximal soil sensors for fine-scale soil attribute mapping in a tropical landscape fosters the research and on-farm operational use of proximal soil sensors in tropical soils.

4. Conclusions

Regression models from a combination of proximal soil sensors produced better soil property predictions than those from individual sensors, excluding the pXRF sensor. In turn, the pXRF sensor produced the best predictions for chemical and physical soil properties among individual sensors and outperformed fusion of the other sensors for soil OC content, clay content, and BD prediction. It is one of the most portable sensors, making it ideal for in situ measurements, and directly measures other soil constituent elements that may be of interest, including plant macro- (K) and micronutrients (Mn, Fe, Zn), heavy metals and pollutants (Cr, Ba, Pb), elements present in minerals linked to soil-forming processes (Ti, Zr), and others. The proximally-sensed eTh content from the gamma-ray sensor showed 0.78 correlation with clay, 0.67 correlation with θ, and moderate correlations (r > 0.50) with OC, SB, and CEC, and is the easiest to operate in the field.

Proximal soil sensor fusion produced better soil property maps through OK of sensor-based predictions compared to OK of the raw data, except for OC. The soil spatial patterns, which in turn respond to relief patterns and processes and water dynamics, were captured with better detail and higher accuracy, by the individual or fused soil sensor data. In situ sensor measurements are fast and do not necessarily require soil sampling and handling (transportation, storage, preparation, laboratory analysis), and thus, can be taken with minimum effort in more sites to increase the sample size and improve spatial coverage and extent. In the new sampling sites, there is no need to take soil samples and run laboratory analysis. On one side, this reduces the time and cost of sampling and data generation, and on the other, it improves the quality of high-resolution soil property maps.