Multi-Sensor Approach for Tropical Soil Fertility Analysis: Comparison of Individual and Combined Performance of VNIR, XRF, and LIBS Spectroscopies

Rapid, cost-effective, and environmentally friendly analysis of key soil fertility attributes requires an ideal combination of sensors. The individual and combined performance of visible and near infrared (VNIR) diffuse reflectance spectroscopy, X-ray fluorescence spectroscopy (XRF), and laser-induced breakdown spectroscopy (LIBS) was assessed for predicting clay, organic matter (OM), cation exchange capacity (CEC), pH, base saturation (V), and extractable (ex-) nutrients in tropical soils. A set of 102 samples, collected from two agricultural fields, with broad ranges of fertility attributes were selected. Two contrasting data fusion approaches have been applied for modeling: (i) merging spectral data of different sensors followed by partial least squares regression (PLS), known as fusion before prediction; and (ii) applying the Granger and Ramanathan (GR) averaging approach, known as fusion after prediction. Results showed VNIR as individual technique to be the best for the prediction of clay and OM content (2.61 ≤ residual prediction deviation (RPD) ≤ 3.37), while the chemical attributes CEC, V, ex-P, ex-K, ex-Ca, and ex-Mg were better predicted (1.82 ≤ RPD ≤ 4.82) by elemental analysis techniques (i.e., XRF and LIBS). Only pH cannot be predicted regardless the technique. The attributes OM, V, and ex-P were best predicted using single-sensor approaches, while the attributes clay, CEC, pH, ex-K, ex-Ca, and ex-Mg were overall best predicted using multi-sensor approaches. Regarding the performance of the multi-sensor approaches, ex-K, ex-Ca, and ex-Mg, were best predicted (RPD of 4.98, 5.30, and 4.11 for ex-K, ex-Ca and ex-Mg, respectively) using two-sensor fusion approach (VNIR + XRF for ex-K and XRF + LIBS for ex-Ca and ex-Mg), while clay, CEC and pH were best predicted (RPD of 4.02, 2.63, and 1.32 for clay, CEC, and pH, respectively) with the three-sensor fusion approach (VNIR + XRF + LIBS). Therefore, the best combination of sensors for predicting key fertility attributes proved to be attribute-specific, which is a drawback of the data fusion approach. The present work is pioneering in highlighting benefits and limitations of the in tandem application of VNIR, XRF, and LIBS spectroscopies for fertility analysis in tropical soils.


Introduction
Spectro-analytical techniques that allow direct analysis of samples are promising for environmentally friendly soil fertility characterization [1]. Integrating these sensing techniques with laboratory analytical methods in hybrid laboratories [2] can modernize the traditional methods for soil analysis [3] that will enable a new paradigm in soil management. The maturation of direct analysis techniques-i.e., those requiring minimal or no sample preparation-would also enable more accurate mapping of agricultural fields by means of in situ or mobile laboratory analysis [4]. Proximal soil sensing technologies are also promising for precision agriculture [5] and pedometrics [6] approaches.
For optimal fertilizer management, spatio-temporal information of key soil fertility attributes as clay, organic matter (OM), cation exchange capacity (CEC), pH, base saturation (V), and extractable (ex-) nutrients is required [4]. Although studies using single sensing techniques have demonstrated the potential to assess some of the above-mentioned attributes [7], it is unlikely that a single sensor can completely characterize the soil fertility to a sufficient degree of accuracy [5,8]. Each sensor presents a perspective regarding the possibilities of predicting soil attributes, which is related with the fundamentals of each technique and its operation mode. Thus, multi-sensor approaches allow the integration of information related to distinct soil properties, yielding more comprehensive and accurate characterizations of key soil fertility attributes [9][10][11].
Visible and near infrared diffuse reflectance spectroscopy (VNIR), X-ray fluorescence spectroscopy (XRF), and laser-induced breakdown spectroscopy (LIBS) are multi-source techniques compatible with direct analysis of soil samples. The visible part of VNIR spectra is sensitive to valence electron excitations of some atoms and functional groups, while the near infrared part responds to non-fundamental vibrations of specific molecules [12]. In soil analysis, the spectral features present in the VNIR region are strictly related to mineralogical and organic components [13]. In turn, elemental analysis techniques, as XRF and LIBS, provide the characterization of the elemental constitution of a wide variety of chemical elements present in soil samples [14]. XRF and LIBS are also considered complementary [14]. For example, elements with atomic number (Z) lower than 12 are hardly detected by XRF technique; however, elements such C (Z = 6) and Na (Z = 11) can be detected by LIBS. Hence, this study hypothesizes that the combination of data provided by VNIR sensors with those from elemental analysis sensors-i.e., XRF and LIBS-might provide complete characterization of key soil fertility attributes [8,15].
Although VNIR, XRF, and LIBS techniques have been evaluated individually for predicting fertility attributes [16][17][18][19][20][21], the combination of these techniques is still at its early stages of development. Recent studies have evaluated different data fusion approaches for combining VNIR and XRF data for predicting fertility attributes [22][23][24]. To the best of our knowledge, only Xu et al. [11] combined VNIR, XRF, and LIBS sensors for predicting soil fertility attributes (OM, pH, ex-P, ex-K, ex-N, and total N). Although these studies show the synergy among these techniques, they do not make a comprehensive evaluation of all the key soil fertility attributes. Moreover, the combination of these three sensors have never been evaluated for the analysis of key fertility attributed in tropical soils. Therefore, the present study aimed at assessing the performance of the individual and combined use of VNIR, XRF, and LIBS sensors for predicting key soil fertility attributes in tropical soils. Two contrasting approaches of data fusion modeling has been applied: (i) a front-end approach that combines the raw data of each sensor, followed by partial least squares (PLS); and (ii) a back-end approach, that applies an adapted approach from the Granger and Ramanathan (GR) averaging method. While the former method can be called fusion before prediction, the latter can be designated as fusion after prediction.

Soil Samples and Reference Analysis
A total of 102 samples were selected from the soil sample bank of the Precision Agriculture Laboratory (LAP) of the University of São Paulo (USP). These samples were selected because their chemical analyses conducted prior to this study yielded broad ranges of variation in the studied soil fertility attributes. These samples were from two agricultural fields with soils characterized as Lixisol and Ferralsol [25], which are common soil types found in Brazil's tropical areas [26]. The coordinates of the Lixisol field and Ferrasol field are 22 • 41 57" S and 47 • 38 33" W and 14 • 06 05" S and 57 • 45 58" W, respectively. The samples were collected from 0-20 cm depth, and stored after air-drying for 48 h and sieving (≤2 mm). The samples were again subjected to laboratory analyses, which provided the results of the reference data (y-variables) that were used to build the predictive models using the different spectral datasets. The contents of clay, OM, CEC, pH, V, ex-P, ex-K, ex-Ca, and ex-Mg were determined following the methods described by Van Raij et al. [27], in which, extractable nutrients were determined using ion exchange resin extraction; OM content by oxidation with potassium dichromate solution, pH by calcium chloride solution, and clay by using the Bouyoucos hydrometer method. The CEC was calculated by the sum of the soil potential acidity determined via buffer solution (SMP) plus the sum of bases. The V was calculated by the ratio between the sum of bases and CEC. Interrelationships between different soil attributes were calculated using Pearson correlation matrix.

Sample Preparation
The data acquisition with XRF and VNIR sensors was performed on the dry and sieved samples. For the XRF measurements, about ten grams of each loose soil sample was placed in a polyethylene cup of 31-mm in diameter (Chemplex Industries Inc., Palm City, FL, USA), whose bottom was sealed with a 4-µm thick polypropylene film (SPEX CertiPrep Inc., Metuchen, NJ, USA). For the LIBS analysis, approximately eight grams of each soil sample was mixed (10% w w −1 ) with a binder material (microcrystalline cellulose, Merck, Darmstadt, Germany), grounded, and homogenized using a planetary ball mill (PM 200 mill, Retsch, Haan, Germany) for 20 min before they were finally pelletized.

VNIR Measurements and Spectra Pre-Processing
The spectrometer of Veris MSP3 (Veris Technologies, Salina, Kansas, USA) was used in benchtop mode to acquire the VNIR spectra. This system uses a tungsten halogen lamp as the energy source and two detection systems, a CCD array (USB4000, Ocean optics, Largo, FL, USA) and an InGaAs photodiode-array (C9914GB, Hamamatsu Photonics, Hamamatsu, Japan), to register spectra from 343 to 2222 nm, with ±5 nm of spectral resolution. The same spectrometer used in this study was described elsewhere by Knadel et al. [28] and Debaene et al. [29]. As a quality control, the instrument automatically checks its spectrum intensity using four factory reference materials during its initialization. The sensor was self-calibrated, by making a dark and white reference measurements, before each spectra acquisition. For laboratory scanning, the samples were placed in a sample holder that prevents ambient light. Each sample was scanned in triplicate, by repositioning the sample after each reading to account for the micro heterogeneity of the sample [30]. Before the spectra pre-processing steps, the replicates were averaged.
The raw spectra were first reduced from 343-2222 nm to 437-2149 nm range due to noise at the spectrum edges. Thereafter, the link between the two detectors, an artifact at 1040 nm, was corrected following the method proposed by Mouazen et al. [31]. These spectra were then pre-processed following the successive order: standard normal variate (SNV) > maximum normalization > first derivative with Savitzky-Golay > smoothing with Savitzky-Golay. These four pre-processing algorithms are frequently adopted in the literature. The SNV is a scattering correction method that is applied to soil VNIR spectra to remove the multiplicative interferences of particle size [32]. The maximum normalization was carried out to get all spectral data to approximately the same scale [33], while the first derivative was applied to improve the signal-to-noise ratio by enhancing weak spectral features [34]. The smoothing approach was used to remove noise and improve the signalto-noise ratio [35,36] further. All data pre-processing steps were performed using the Unscrambler ® version 10.5.1 (Camo AS, Oslo, Norway).

XRF Measurements and Spectra Pre-Processing
A portable device Tracer III-SD model (Bruker AXS, Madison, Wisconsin, EUA) was used for XRF data acquisition. This instrument has a 4 W Rh X-ray tube and a Peltier-cooled Silicon Drift Detector with 2048 channels. The X-ray tube was configured for voltage and current of 35 kV and 7 µA, respectively [18]. No filter was used and the scanning time was 90 s, which was performed under atmospheric pressure [18]. Three readings were taken from each soil specimen at three different spots, and these were then averaged to be considered for spectra pre-processing and modeling.
The acquired spectra were normalized by the detector live time and evaluated in counts of photons per second (cps). Considering the area under each peak, 12 spectral lines (K-lines of Al, Si, K, Ca, Ti, Mn, Fe, Ni, and Cu, and the scattering peaks Rh-Lα Thomson, Rh-Kα Compton, and Rh-Kα Thomson) were selected to be used as explanatory variables, following the criteria recommended by Tavares et al. [37]. Finally, the K-lines of Al, Si, K, Ca, Ti, Mn, Fe, Ni, and Cu were normalized by the Compton peak. XRF emission lines were obtained using Artax ® software (Bruker AXS, Madison, USA), whose pre-processing was conducted in Excel 2016 (Microsoft Corporation, Redmond, WA, USA).

LIBS Measurements and Spectra Pre-Processing
A benchtop LIBS system composed by a pulsed Nd:YAG laser at 1064 nm, generating 5 ns pulses of up to 365 mJ (Brilliant, Quantel, France) and an ESA 3000 spectrometer (LLA Instruments GmbH, Berlin, Germany) was used. The laser was focused on the sample surface by a convergent lens. Pressed pellets were placed into a plastic sample holder positioned in a two axes manually-controlled translation stage, movable in the plane orthogonal to the laser direction. A laminar stream of argon (5.0 L min −1 ) was continuously fed from the bottom of the sample holder in order to dislocate the atmospheric air around the sample surface. The emission from the plasma was collected by using a telescope composed of 50 mm and 80 mm focal length fused silica lenses and coupled to the entrance slit of the spectrometer using an optical fiber. The collection angle with respect to the laser optical axis was 25 • . This spectrometer registers spectra from 200 to 780 nm with a resolution oscillating from 5 pm at 200 nm to 19 pm at 780 nm (totaling 53,717 variables). The same LIBS system used in this study was described elsewhere by Nunes et al. [38].
LIBS experimental conditions were carried out with 65 mJ laser pulses, 19.5 cm of lensto-sample distance (180 µm spot size providing 255 J cm −2 laser fluence), 15 accumulated laser pulses, 2 µs of delay time, and 7 µs of integration time gate. These instrumental conditions were optimized in initial tests to obtain the maximum signal-to-noise ratio of the emission lines of interest, as suggested by Nunes et al. [39]. In order to consider the analytes micro-heterogeneity in the test samples, 21 laser shots were applied at different positions of the pellet surface, and the corresponding analytical signals were averaged in one reading to be used for the following pre-processing and modeling steps.
Initially, the LIBS spectra were pre-processed by downscaling them from 200 to 540 nm, which is a region containing the most valuable and high-intensity emission lines [40]. Afterwards, the interval successive projection algorithm (iSPA) was applied to select significant variables [41]. To employ this algorithm, the LIBS spectra from 200 to 540 nm were divided into 160 intervals of 243 variables. The LIBS spectra were pre-processed in MATLAB (Math-Works Inc., Natick, MA, USA) using the script developed by Gomes et al. [41]. The regions selected for each fertility attribute are shown in Table A1 (Appendix A). The mean raw spectra acquired with the VNIR, XRF, and LIBS technique are shown in Figure 1. Section). The mean raw spectra acquired with the VNIR, XRF, and LIBS technique are shown in Figure 1.

Modeling
After pre-processing, spectral variables were associated with the laboratory measured soil attributes by establishing predictive models. This was done using the spectral data of each single-sensor and for all possible combinations of sensors, through the following associations: VNIR + XRF, VNIR + LIBS, XRF + LIBS, and VNIR + XRF+ LIBS. The calibration and validation subsets comprised, respectively, of 70% and 30% of the total number of samples, which were divided using the Kennard-Stone algorithm [42] executed on the measured soil attributes. The prediction performance was evaluated through the determination coefficient (R 2 ), root mean square error (RMSE), and the residual prediction deviation (RPD). Based on the RPD values, the prediction quality of developed models were classified into four classes adapted from Chang et al. [43]: poor models (RPD < 1.40), reasonable models (1.40 ≤ RPD < 2.00), good models (2.00 ≤ RPD < 3.00), and excellent models (RPD ≥ 3.00). The relative improvement (RI) of the predictions accomplished by means of the multi-sensor approaches were also calculated, in terms of percentage of RMSE [44,45], and compared to the best prediction obtained using the single-sensor (VNIR, XRF, and LIBS, individually) and two-sensor approach (VNIR + XRF, VNIR + LIBS, and XRF + LIBS).

Single-Sensor Predictive Models
The predictive models using VNIR and LIBS data individually were established by means of PLS regression [16,21]. The number of latent variables adopted for each PLS model was determined for the model in cross-validation that resulted in the lowest RMSE. For the XRF, models were calibrated with multiple linear regression (MLR), as suggested by Tavares et al. [37]. All the calibrations and validations were performed using the Unscrambler software, version 10.5.1 (Camo AS, Oslo, Norway).

Data Fusion Approaches
Two data fusion approaches were used to build the predictive models using multisensor data. The first was designated as spectral fusion (SF) method. It is a front-end approach, consisting of merging the pre-processed spectral data of each combination of sensor (VNIR + XRF, VNIR + LIBS, XRF + LIBS, and VNIR + XRF + LIBS) in one matrix. The merged data were scaled by the standard deviation to create an even distribution of variances, before they were subjected to PLS regression to develop a model for each soil attribute. The number of latent variables for each PLS model was determined according to the leave-one-out cross-validation that resulted in the lowest RMSE.
The second data fusion method was adapted from the Granger and Ramanathan (GR) averaging method [46], as performed by Tavares et al. [9]. In the present study, the GR approach consisted of calibrating a multiple linear regression (MLR) with independent variables being the predictions made with the sensors, plus the prediction made with the multi-sensor merged data (i.e., SF approach). For example, for the combination between XRF and VNIR sensors, the GR method was applied according to Equation (1) where Y is the soil property to be estimated; X VNIR and X XRF are the predictions made by VNIR and XRF single-models, respectively; X SF is the prediction made by the SF approach using the VNIR and XRF merged data; W 0 , W VNIR , W XRF , and W SF are the parameters of the MLR determined by least squares method, where the first parameter is the value of the line intercept and the others are the weights of the predictions.
Equation (1) was also applied for combinations of LIBS + VNIR, LIBS + XRF, and VNIR + XRF + LIBS. All data fusion models were calibrated and validated using the Unscrambler software, version 10.5.1 (Camo AS, Oslo, Norway).

Characterisation of the Laboratory Measured Soil Properties
The boxplots in Figure 2 show that the range and standard deviation (SD) of soil attributes in the calibration set are close to those in the validation set, which is expected since we applied Kennard Stone method in sample split to avoid undesirable influences on the prediction accuracy that are not sensor-related [12,47]. If different range and SD exist between the validation and the calibration sets, this would result in poor prediction of concentrations that are not accounted for in the calibration set.

Characterisation of the Laboratory Measured Soil Properties
The boxplots in Figure 2 show that the range and standard deviation (SD) of soil attributes in the calibration set are close to those in the validation set, which is expected since we applied Kennard Stone method in sample split to avoid undesirable influences on the prediction accuracy that are not sensor-related [12,47]. If different range and SD exist between the validation and the calibration sets, this would result in poor prediction of concentrations that are not accounted for in the calibration set.  Figure 3 shows the correlation matrix of the soil fertility attributes. These interrelationships aid to understand reasons of successful determination of soil attributes having only indirect spectral signatures in the studied sensors. For example, pH or V has no emission lines in XRF and LIBS spectra, nor direct spectral signature in the VNIR range. In this study, it is important to highlight the following interrelationships: (i) CEC, V, ex-Ca, and ex-Mg were closely related with one another, showing Pearson correlations ranging from 0.85 to 0.94; (ii) ex-K had strong correlations with clay and V (r = 0.73 for both); (iii) clay showed strong correlations with V, ex-K, and ex-Ca (0.73 ≤ r ≤ 0.82); (iv) OM content was moderately correlated with all attributes (0.44 ≤ r ≤ 0.62), except for ex-P and pH (r of 0.06 and −0.02, respectively); (v) pH presented moderate correlations with CEC, V, ex-K, ex-Ca, and ex-Mg (0.40 ≤ r ≤ 0.50); and (vi) ex-P was the attribute that presented weaker correlations (−0.23 ≤ r ≤ 0.06), indicating P as an independent attribute from the other attributes, with potential of poor prediction using spectroscopic methods that do not present spectral response for P (e.g., VNIR and XRF).  Figure 3 shows the correlation matrix of the soil fertility attributes. These interrelationships aid to understand reasons of successful determination of soil attributes having only indirect spectral signatures in the studied sensors. For example, pH or V has no emission lines in XRF and LIBS spectra, nor direct spectral signature in the VNIR range. In this study, it is important to highlight the following interrelationships: (i) CEC, V, ex-Ca, and ex-Mg were closely related with one another, showing Pearson correlations ranging from 0.85 to 0.94; (ii) ex-K had strong correlations with clay and V (r = 0.73 for both); (iii) clay showed strong correlations with V, ex-K, and ex-Ca (0.73 ≤ r ≤ 0.82); (iv) OM content was moderately correlated with all attributes (0.44 ≤ r ≤ 0.62), except for ex-P and pH (r of 0.06 and −0.02, respectively); (v) pH presented moderate correlations with CEC, V, ex-K, ex-Ca, and ex-Mg (0.40 ≤ r ≤ 0.50); and (vi) ex-P was the attribute that presented weaker correlations (−0.23 ≤ r ≤ 0.06), indicating P as an independent attribute from the other attributes, with potential of poor prediction using spectroscopic methods that do not present spectral response for P (e.g., VNIR and XRF).

Prediction Performances of Single-Sensor and Data Fusion Approaches
The prediction results of individual and combined application of VNIR, XRF, and LIBS sensors are presented in Table 1 (see details in Table A2, in the Appendix Section). Considering the single-sensor prediction performance, VNIR presented excellent and good prediction performance for clay (RPD = 3.37), OM (RPD = 2.61), and V (RPD = 2.26). Among the three sensors, VNIR showed the lowest RMSE for clay and OM prediction (Table A2). The XRF sensor showed good prediction performance for clay (RPD = 3.  (Table A2). The pH was the only attribute that could not be satisfactorily predicted (RPD < 1.12) by any of the single-sensor approaches, while the ex-P was only satisfactorily predicted by LIBS (RPD = 1.82).

Prediction Performances of Single-Sensor and Data Fusion Approaches
The prediction results of individual and combined application of VNIR, XRF, and LIBS sensors are presented in Table 1 (see details in Table A2, in the Appendix A). Considering the single-sensor prediction performance, VNIR presented excellent and good prediction performance for clay (RPD = 3.37), OM (RPD = 2.61), and V (RPD = 2.26). Among the three sensors, VNIR showed the lowest RMSE for clay and OM prediction (Table A2) (Table A2). The pH was the only attribute that could not be satisfactorily predicted (RPD < 1.12) by any of the single-sensor approaches, while the ex-P was only satisfactorily predicted by LIBS (RPD = 1.82). Table 2 shows the RI for the predictions performed using multi-sensor data through the two applied data fusion approaches (SF and GR), compared with the best single-sensor approach. This indicator shows the improvement (when positive) or the deterioration (when negative) obtained in the sensor fusion relative improvement [44,45]. The attributes OM, V, and ex-P were better predicted by the individual than combined sensors. All the other attributes (clay, CEC, pH, ex-K, ex-Ca, and ex-Mg) were better predicted by multi-sensor approaches. Table 1. Residual prediction deviation (RPD) obtained for the validation set (n = 34) of clay, organic matter (OM), cation exchange capacity (CEC), pH, base saturation (V), and extractable (ex-) nutrients (ex-P, ex-K, ex-Ca, and ex-Mg) using the visible and near infrared (VNIR), X-ray florescence (XRF), and laser-induced breakdown spectroscopy (LIBS) data individually and different combination of them through spectra fusion (SF) and Granger  RPD values for the same soil attribute were compared and presented on green scale, highlighting the highest values within each soil attribute. The root-mean-square error (RMSE) are given in g dm −3 for clay and OM; in mmol c dm −3 for CEC, ex-K, ex-Ca, and ex-Mg; in % for V; and, and in mg dm −3 for ex-P. Results of coefficient of determination (R 2 ) and RMSE are detailed in the Appendix A (Table A2). Table 2. Relative improvement (RI) and root-mean-square error (RMSE) achieved for each multi-sensor approach using spectra fusion (SF) and Granger and Ramanathan (GR) data fusion methods, in comparison with the best prediction obtained using the single-sensor approaches. The RI values for clay prediction ranged from −11.1% to 16.3%, for the SF approach, and from −3.2% to 15.2%, for the GR approach. For CEC prediction, the RI ranged from −31.1% to 1.4%, using SF, and from −32.0% to 2.3%, using GR. Although the pH models did not perform satisfactorily, their RI values ranged from −6.5% to 14.8% and 1.9% to 15.3% for SF and GR, respectively. For the ex-K prediction, the RI ranged from −117.4% to −12.6% for SF, and −104.5% to 14.4% for GR. For ex-Ca, the RI varied between −70.2% and 6.2% for SF, and between −23.9% and 9.2% for GR. For ex-Mg prediction, the RI values ranged from −28.5% to 24.2%, using SF, and from −12.4% to 27.2%, using GR. VNIR + XRF led to positive improvement in the prediction of clay, pH and ex-K (with GR only), although improvements were larger in VNIR + LIBS for clay and pH only. As expected, the fusion of XRF and LIBS improved the prediction of the extractable nutrients mainly, ex-K, ex-Ca, and ex-Mg, with minor improvement in CEC, with SF, and pH, with GR ( Table 2). The fusion of the three sensors have scored the largest number of improvements with positive RI values for all studied attributes except for OM, V, and ex-P. In general, the performance of the back-end data fusion strategy using the GR approach outperformed the SF approach in 72.2% of cases. In summary, the improvement in predictive performance varied according to the combination of sensors as well as the data fusion strategy used; there was no one optimal approach for the prediction of all the attributes analyzed. Table 3 compares the prediction performances of the three-sensor approach with the best performances obtained with the two-sensor and the single-sensor approaches. Results revealed that the best predictions of OM, V, and ex-P were obtained using the single-sensor approach (with VNIR for OM, and with LIBS for V and ex-P). Predictions of these attributes using two sensors showed degradation in predictive performance, with the best two-sensor approach showing an RMSE increment of 7.1%, 0.7%, and 18.2% for OM, V, and ex-P, respectively. A similar deterioration in the prediction performance of these attributes (OM, V, and ex-P) was also observed by the three-sensor approach, with an increase in the prediction error oscillating between 6.9% and 74.2% in comparison to the best performing single-sensor approach.

Clay
Concerning ex-K, ex-Ca, and ex-Mg, better prediction results were obtained when using the two-sensor approach. The VNIR + XRF combinations were the most promising for ex-K, while XRF + LIBS stood out for ex-Ca and ex-Mg. The performance increment obtained by the above-mentioned two-sensor approaches, relatively to the best singlesensor strategy, was 14.4%, 9.2%, and 27.2% for ex-K, ex-Ca, and ex-Mg, respectively. When combining a third sensor, the predictive performance was degraded, with a reduction of 2.9%, 6.5%, and 6.8% for ex-K, ex-Ca, and ex-Mg, respectively, relative to the best two-sensor approach.
Concerning the clay, CEC, and pH the best prediction performances were recorded by the three-sensor approach, which reduced the error to a range between 2.3% and 16.3%, in comparison to the best single-sensor approach, and between 0.9% and 2.2%, relative to the best approach using two sensors.
Due to the presence of synergistic information in the data, improved predictions were obtained for a large number of attributes (six out of nine) when using the three-sensor approach. The combination between XRF and LIBS showed synergy for the prediction of five out of nine attributes, while both the combinations VNIR + LIBS and VNIR + XRF showed synergy for three out of nine attributes. Table 3. Relative improvement (RI) achieved for the predictions using the best three-sensor approach in comparison with the best prediction obtained using the two-sensor, and the single-sensor approaches. The root mean square error (RMSE) of each approach was also presented.  1 Organic matter; 2 cation exchange capacity; 3 base saturation; 4 extractable (ex-) nutrients; 5 related technique or combination of techniques; 6 related data fusion approach; 7 RI compared with the best single-sensor approach; 8 RI compared with the best two-sensor approach; 9 visible and near infrared diffuse reflectance spectroscopy; 10 X-ray fluorescence spectroscopy; 11 laser-induced breakdown spectroscopy; 12 Granger and Ramanathan approach; and 13 spectral fusion approach. For each soil attribute, bolded RMSE values indicate the approach with the lowest prediction error.

Discussion
Soil spectral analyses via the tandem use of VNIR, XRF, and LIBS sensors have two potential benefits: (i) extend the set of predicted soil attributes in comparison with single-sensor approaches; and (ii) improve the models' predictive performance by exploring synergy in the data through data fusion approaches, both are discussed in Sections 4.1 and 4.2, respectively.

Individual Performance of VNIR, XRF, and LIBS Sensors: Coverage of Key Soil Fertility Attributes
The results observed in this work show that the tandem application of VNIR, XRF, and LIBS techniques, allows extending the number of predicted fertility attributes with optimal performances (Table A3, in the Appendix A). For individual sensor, limited number of attributes could be detected and optimal predictions are technology specific. For example, excellent and good predictions could be achieved for clay and OM, respectively, with both the VNIR and LIBS sensors; although, the VNIR sensor allowed more accurate predictions of these attributes, since both attributes have direct spectral responses in the NIR range [48]. A similar trend was observed for the chemical attributes CEC, V, ex-K, ex-Ca, and ex-Mg, whose optimal prediction performances were obtained with elemental analysis sensors (with RPD oscillating between 2.57 and 4.82) due to the fact these attributes are closely related to the emission lines present in these sensors.
The LIBS spectra presented emission lines for C, P, and Mg (C I 247.856 nm, P I 213.6182 and P I 214.914 nm, Mg I 277.983, and Mg I 285.213 nm), elements not observable in the corresponding XRF spectra. The presence of these emission lines explains the better performance of LIBS for OM, ex-P, and ex-Mg prediction than XRF. The XRF spectra showed K emission line (in 3.31 keV), which was not observed in the LIBS data, explaining its superior performance for ex-K prediction. The excellent ex-Ca prediction of both elemental techniques is well supported by appropriate detection of Ca emission lines (Ca Kα 3.69 keV for XRF, and Ca II 315.887 nm and Ca II 317.933 nm for LIBS). The presence of Ca emission lines also explains the similar performance for V and CEC prediction obtained by both XRF and LIBS (good predictions for CEC and excellent for V), since both attributes showed strong correlation with the ex-Ca (r ≥ 0.92; Figure 3).
The attributes with the worst predictive performance were the ex-P and the pH. Regarding ex-P, LIBS was the only sensor to present a satisfactory prediction for this attribute. Although LIBS shows well defined emission lines for P, total contents of P may not always correlate with ex-P contents due to the complex chemistry of this element in tropical soils [49]. P in tropical soils is commonly associated, in a non-extractable manner, with Fe and Al oxides [50], which makes the modeling of ex-P via elemental analysis techniques challenging [18,19]. A possible alternative to improve the predictive performance of ex-P and pH models is through the integration of electrochemical sensors (e.g., ion-selective electrodes (ISE) and the ion-sensitive field-effect transistor (ISFET)) to the multi-sensor array used, since with the former technology ensure direct soil analysis of the extractable contents of some ions (e.g., H + and PO 4 3− ) [4,7]. In general, the predictive capability of key soil fertility attributes decreased as follows: LIBS > XRF > VNIR. These results differed from those presented by Xu et al. [11], which is the only research that explored VNIR, XRF, and LIBS sensors for predicting soil fertility attributes. The authors obtained better prediction performances with the VNIR sensor, in comparison to the XRF and LIBS sensors, for OM, pH, ex-P, ex-K, and ex-N. Our results, however, are in line with those of other research carried out in tropical soils that also evaluated the standalone use of VNIR [2,[51][52][53][54] and XRF [19,[55][56][57][58] techniques. Studies using LIBS sensors for the assessment of fertility attributes are incipient in tropical soils. However, the satisfactory performance obtained for clay is supported by that obtained by Villas-Boas et al. [59]. In turn, our poor pH predictions contrasted with the satisfactory predictions obtained by Ferreira et al. [60].

Combined Performance of VNIR, XRF, and LIBS Sensors: Synergy for Predicting Key Soil Fertility Attributes
The prediction of OM, V, and ex-P did not show any synergy by the multi-sensor approaches evaluated in this study, as they were better predicted by each sensor individually. However, improved prediction performance was obtained for all the other attributes (clay, CEC, pH, ex-K, ex-Ca, and ex-Mg), with RI ranging between 0.2% and 27.2% when using at least one of the multi-sensor approaches. The absence of synergy in the information provided by VNIR and elemental analysis sensors has already been reported in the literature for pH, ex-K, and OM [9,11]. Conversely, O'Rourke et al. [44], assessing the combination XRF + VNIR in Australian soils, reported RI ranging from 15% to 44% for clay, CEC, pH, ex-K, ex-Ca, and ex-Mg. Also evaluating data fusion between VNIR and XRF sensors, Zhang and Hartemink [24] obtained RI values oscillating from 3% to 20% for the prediction of clay, pH, and total carbon in North American soils. A similar study achieved RI of 26% for CEC in Chinese soils [22]. Also in Chinese soils, Xu et al. [10] reached an RI varying between 4% and 35% for OM prediction by the combination of LIBS and attenuated total reflectance Fourier-transform mid-infrared spectroscopy (FTIR-ATR).
In general, the prediction quality of the studied attributes by the sensor fusion decreased in the order of: VNIR + XRF + LIBS > XRF + LIBS > VNIR + LIBS > VNIR + XRF. Our results showed that there was no unique optimal sensor combination for predicting all the key soil fertility attributes, and that this is attribute specific. A similar behavior was also observed by Xu et al. [11]. Thus, the performance of different combinations of sensors should be evaluated by trial-and-error and that it is too soon to establish a clear trend, which requires further research. This means that data fusion is in its infancy stage, and more works will be necessary to recommend the optimal set of sensors for one or a set of soil attribute.
This study was conducted using 102 soil samples with broad variation of fertility attributes. They belong to Lixisols and Ferralsols, which are representative and common types of soil in Brazilian tropical agricultural lands. In addition, this pioneering evaluation provided useful information to help PSS users to understand the advantages and drawbacks of the combined use of VNIR, XRF, and LIBS sensors for soil fertility analysis. Nevertheless, future studies should consider larger datasets that account for other types of soils, soil mineralogy, textural classes, concentration range, and different agricultural practices.

Operational Aspects Related to Sample Preparation
It is important to consider practical aspects of sample preparation necessary to ensure the application of each technique. The VNIR and XRF techniques are less demanding in terms of sample preparation, and can be applied on loose soil samples (i.e., dried and sieved with 2 mm of particle size) without major impacts on the sensors performance [30,51].
However, the LIBS demands further sample preparation steps including particle size reduction with grinding, the addition of binding agent, followed by pressing to create pellets that are necessary for optimal LIBS performance. Pelletizing promotes a more reproducible laser-sample interaction that is necessary to maintain the stoichiometric proportion of the craters formed by the laser pulses [61]. These extra sample preparation steps to ensure the collection of reliable data make it challenging to apply LIBS to run analyses directly in the field [62]. However, despite the extra time consigned to sample preparation, it would not make it unfeasible to use LIBS in controlled environments (e.g., hybrid laboratories and mobile laboratories) that enable the execution of such sample conditioning [4].

Conclusions
This study shows that the fusion of data from visible and near infrared diffuse reflectance spectroscopy (VNIR), X-ray fluorescence spectroscopy (XRF), and laser-induced breakdown spectroscopy (LIBS) can extend the number of fertility attributes predicted with optimal performance, as well as improve the prediction accuracy by exploiting the synergy through data fusion techniques.
Comparing among the individual sensor performance, the VNIR technique performed best for the prediction of clay and organic matter (OM) (2.61 ≤ RPD ≤ 3.37), while the elemental analysis techniques showed better performance for the prediction of cation exchange capacity (CEC), base saturation (V), and the extractable (ex-) K, Ca, and Mg (2.61 ≤ RPD ≤ 3.37). The LIBS stood out for V, ex-P and ex-Mg prediction, while the XRF resulted in optimal predictions for CEC, ex-K, and ex-Ca. The attributes with the worst predictive performance were the ex-P (0.80 ≤ RPD ≤ 1.82) and pH (RPD ≤ 1.12), with the latter performing poorly in all individual and combined sensor options tested in this study.
Regarding the performance of the multi-sensor approaches, ex-K, ex-Ca, and ex-Mg, were better predicted when using two sensors (VNIR + XRF for ex-K, and XRF + LIBS for ex-Ca and ex-Mg), while clay, CEC and pH were best predicted with the three-sensor (VNIR + XRF + LIBS). However, the combined use of sensors did not always lead to improvement in the prediction results. For instance, the best sensing method for OM, V, and ex-P prediction was the single-sensor approach, showing that there was no synergy between the evaluated sensors for the prediction of these attributes. Furthermore, the best combination of sensors for predicting key fertility attributes proved to be attribute-specific, which is a drawback of the in tandem application of these analytical techniques. Another drawback to integrate the addressed techniques is the need for different sample preparations, where LIBS needs extra steps to conform the loose soil samples into pellets.
Finally, the present work evidenced benefits and limitations of the in tandem application of VNIR, XRF, and LIBS spectroscopies for fertility analysis in tropical soils. Further research should be encouraged to improve analytical protocols using multi-sensor approaches, to enable practical, environmentally friendly, and accurate analysis, allowing its incorporation into future hybrid laboratories.

Acknowledgments:
We would like to thank the technician Eduardo de Almeida, from the Center for Nuclear Energy in Agriculture of the University of São Paulo (USP-CENA), for the support with X-ray fluorescence equipment; and Fábio L. Melquiades and Felipe R. dos Santos, both from the Londrina State University (UEL), for the support with the LIBS data pre-processing.

Conflicts of Interest:
The authors declare no conflict of interest.
Appendix A Table A1. Intervals of variables selected by the interval successive projection algorithm (iSPA) in the LIBS spectra for the predictive modeling of each soil fertility attribute.  Table A2. Prediction results of the validation set (n = 34) obtained using single VNIR, XRF, and LIBS data alone and using multi-sensor data, combined through spectra fusion (SF) and Granger and Ramanathan (GR) approaches.  1 Organic matter; 2 cation exchange capacity; 3 base saturation; 4 extractable (ex-) nutrients (ex-P, ex-K, ex-Ca, and ex-Mg). The coefficient of determination (R 2 ) values for the same soil attribute were compared and presented on green scale, highlighting the highest values within each soil attribute. The root-mean-square error (RMSE) are given in g dm −3 for clay and OM; in mmol c dm −3 for CEC, ex-K, ex-Ca, and ex-Mg; in % for V; and, for ex-P, the RMSE was given in mg dm −3 . Table A3. Qualitative interpretation * of the residual prediction deviation (RPD) performance obtained for each of the evaluated approaches using single-sensor (VNIR, XRF, and LIBS) and multi-sensor (VNIR + XRF, VNIR + LIBS, XRF + LIBS, and VNIR + XRF + LIBS).