Multi-Sensor Soil Probe and Machine Learning Modeling

1. Introduction

The need for cost-effective, rapid, deep, and comprehensive soil health characterization in support of climate-smart agricultural management, soil carbon accounting, precision agriculture applications, and digital twins in smart farming is profound [1,2,3]. Traditional approaches to characterize soils are laborious entailing the ex situ collection of soil samples in individual horizons/layers or soil coring, soil analytics in the laboratory, and digital soil mapping and modeling. In situ proximal sensing in the near surface dates to the late 1990’s and early 2000 with Ben-Dor et al. [4] reporting the first instrument of its nature being a soil penetrometer [5] that later was coupled to a window regulating mechanism that collected reflected light, enabling one to view the color and structure of the soil profile [6]. Poggio et al. [7] conducted a laboratory-based evaluation of the optical performance of a soil penetrometer that include visible and near-infrared (VisNIR) optics that acknowledged the contributions of Rooney to the design. Recently, the advent of proximal soil sensor technology and artificial intelligence (AI) soil predictive modeling has excelled to quantify soil health properties, especially soil organic carbon [8]. Soil measurements made in the laboratory under controlled conditions are still considered the “gold” standard in terms of accuracy and precision of measurements, with soil sensors calibrated and validated against these standards. However, the mismatch of the spatial scale and sample support of laboratory and in situ field sensors is stark and has not been sufficiently addressed in research investigations. Sample sizes required to support conventional laboratory analysis often substantially exceed the spatial scale within which soil properties vary. Proximal sensors are in close proximity to soil samples with the potential to continuously characterize the variability of soils along soil profiles [9], while conventional extracted soil samples used for laboratory-based soil analytics are low in volume, mass, and vertical resolution. Figure 1 demonstrates the relationship between in situ proximal sensing and soil sampling for the measurement of soil potassium. The proximal sensor delineates the profile at a sub-centimeter scale which represents how potassium is distributed in nature. However, the laboratory requires 500g of soil sample which equates to about a 30 cm section of a standard core to run a full laboratory testing panel. Of this 500g approximately 50-75g is utilized for nutrient testing. The result from the lab shows one value (55 ppm) while the range detected by the proximal sensor varies from 10 ppm to 98 ppm. Which is more likely to be a true representation of how potassium is distributed in the soil profile? How to compare the two? Which is the “gold” standard?

Another shortcoming of ex situ soil sampling is that the in situ co-relation between soil properties and attributes is disrupted during extraction which further degrades the utility of soil information derived from conventional methods and obfuscates the intent of the survey. For example, when mapping soil-water properties for irrigation management at field-scale, understanding the in situ relationship between grain size and packing, structure, density, microbial gums, and the depth and thickness of soil horizons is critical. By removing a soil sample and performing a laboratory soil texture test, all the corresponding contextual soil information is detached, thus degrading the value of the soil texture data as a predictor of water holding capacity. This is typical of each laboratory testing procedure, whether physical, chemical, or biological. A precise holistic understanding of soil-crop relationships is best informed by in situ soil testing conducted where the roots interact with the soil properties.

In addition to the issues with vertical resolution and loss of sample context, conventional methods are also at a disadvantage with respect to lower spatial resolution. In practice, where what is learned during research is intended to be put to practical use in agriculture, it is not realistic to obtain and test enough ex situ soil samples to create an accurate map with the spatial resolution needed to operate currently available variable rate nutrient and soil amendment applicators and irrigation technologies. The result is that ex situ soil information is interpolated and extrapolated and then digitized in ways that are not recognized or challenged by end users. The resulting ex situ soil data is highly subjective and lacks spatial and information resolution and is not suitable for advanced analytics enabled by a digital twin.

Some understudied research questions are: 1) whether the approach to ground truth in situ field soil sensor data with manually extracted coarser scale ex situ soil samples undergoing laboratory soil analysis realizes the full potential of in situ proximal sensing of soil properties? and 2) which scaling function performs best to link high resolution in situ soil sensor data and coarser resolution ex situ laboratory analytic data?

The most widely studied proximal soil sensors are visible-near-infrared (VNIR) and mid-infrared (MIR) spectral instruments, which have been used to develop soil spectral libraries at global scale [10,11], regional scale [12,13], and national scales; for example, in the U.S. [14,15], Brazil [16], China [17], and Switzerland [18]. The ability to predict soil organic carbon (SOC) using machine learning (ML) with large-scale spectral libraries in the U.S. has shown excellent performance using independent validation data. For example, SOC predictions from VNIR spectra and random forest (RF) modeling achieved a Coefficient of Determination (R2) of 0.95, Ratio of Performance to Inter-Quartile (RPIQ) of 0.81 [14] and R2 of 0.96, RPIQ: 5.18 [19] using rigorous validation assessment for soils in the conterminous U.S. Similarly compelling results in validation mode were achieved for modeling SOC in the U.S. using VNIR spectra and Convolutional Neural Networks (CNN-1) with R2 of 0.83, RPIQ 0.81 and even better results using MIR spectra and CNN-1 R2 of 0.98 and RPIQ of 2.37 [15]. Other physical and chemical soil properties, such as macro- and micro-nutrients, soil texture, cation exchange capacity (CEC), and pH have been predicted widely from diffuse reflectance spectral data [14,20,21,22,23]. In particular, MIR spectral data have fingerprinting capabilities for soil characteristics and elemental content, while VNIR relies on overtones of chemical bonds in spectra (e.g., C-O, C-H, N-H, O-H) [12].

Bulk density (BD) cannot be directly inferred from spectral reflectance data because it relies on associations with other soil properties such as soil texture and SOC. For example, BD (measured using clod-only, core-only, and combined clod and core methods) was predicted using Partial Least Square Regression (PLSR), Cubist, memory-based learner (MBL), and RF from MIR data with R2 in validation mode ranging between 0.64 (PLSR) to 0.81 (MBL) [24]. Davari et al. [25] found that both soil BD (R2 = 0.35) and soil porosity (R2 = 0.16) were poorly predicted using only VNIR spectra suggesting that other sensors such as penetrometers that measure tip and sleeve stress are needed to improve inference capabilities [26,27]. The Soil Condition Analysis System (SCANS) integrates an ex situ soil core scanning system with multiple sensors including a γ-ray attenuation densitometer to measure BD, digital cameras for soil imaging, and VNIR spectrometer [28].

The advantages of spectral soil prediction modeling include that VNIR provides high sample throughput through rapid scanning of samples compared to conventional soil analytics [29,30]. Hyperspectral soil data show significantly higher information content than traditional laboratory soil analytics. Proximal soil sensing is non-destructive and produces no hazardous materials. Another advantage is that once large spectral libraries have been built, they can be reused and improved (e.g., applying novel ML algorithms) over time until they reach model saturation. Review articles of proximal soil sensing technology unequivocally converge in view that proximal soil health sensing of individual soil samples is a mature analytical technique if performed under controlled laboratory conditions using sieving, grinding, and drying (MIR) and sieving and drying (VNIR) of soil samples [31,32,33,34]. Sieving and drying operations are employed to produce comparability among laboratory scanned spectra because soil reflectance spectra are also affected by particle size [35,36,37,38,39] and surface roughness [40,41,42], both of which relate to soil texture.

The emergence of field-based soil spectroscopy using portable or mounted instruments has marked a shift from laboratory settings to in situ field sensing [43]. Some field studies showed significant differences between controlled laboratory and field based VNIR applications due to spatially variable environmental conditions. For example, the study by Hedley et al. [44] used a portable spectroradiometer to predict topsoil SOC from field-moist spectra with a low R2 = 0.39 and Ratio Performance Deviation (RPD) = 1.28 compared to air-dry spectra with R2 = 0.80 and RPD = 2.25 showed significant differences due to the effects of soil moisture. The effects of soil moisture on soil spectral modeling have long been known in the spectral soil community [41,45,46]. According to Seidel et al.'s empirical data [47] (2022), soil moisture effects are more significant in MIR than VNIR applications. Methods, such as external parameter orthogonalization (EPO), direct standardization (DS), global moisture modelling (GMM), slope bias correction (SB), and selective wavelength modelling (SWM), have been suggested to address the application of VNIR under field conditions with varying soil moisture content [19]. In their study dry samples were rewetted with different soil moisture content demonstrating that EPO, DS and GMM account satisfactorily for the effect of moisture in soil spectra. These three methods improved the prediction of SOC substantially with an increase in R2 from almost zero for no correction to over 0.5 and an RPIQ from 0.38 to over 1.7. These findings suggest that the effect of moisture on VNIR modeling of SOC and other soil properties is removable through post-process corrections applied to spectral data. Knadel et al. [48] provided a comprehensive review of mathematical techniques to remove moisture effects from VNIR spectra. However, such approaches are computationally expensive if applied to spectral field data. Data-driven ML methods offer alternatives to removal of soil moisture effects from spectral data by explicitly incorporating moisture data along with spectral and/or other sensor data into soil prediction models.

One such study was presented by Zhou et al. [23] (2024) who analyzed loess soil samples to investigate how changes in soil moisture content impact predictions from VNIR spectra. Various supervised learning and latent variable methods (PLSR; RF; and Support Vector Machines,) were tested with the first derivative-Genetic algorithm (GA)-RF method demonstrating successful predictions of soil moisture with R2 of 0.99 and Relative Prediction Deviation (RPD) of 16.2. Similarly, Lobell & Asner [49] study quantified the strong influence of moisture on spectral reflectance and absorption features. Tan et al. [50] critiqued that many studies using soil spectroscopy focused on dried soil samples in the laboratory under controlled conditions, while techniques to remove soil moisture effects from VNIR spectra are time-consuming and counter-productive in the field. In Tan et al.’s [50] empirical study, soil moisture effects were successfully eliminated from VNIR spectra to model soil organic matter (SOM) using Principal Component Analysis (PCA)-RF coupled with the continuous wavelet transform (CWT). They found that wavelengths at about 580nm, 820nm, and especially the narrow region around 1,400nm are highly correlated regions to SOM using wet soil samples. Validation results to predict SOM from wet samples based on PCA-RF (R2 = 0.84, RPD = 2.53) and dry samples (R2 = 0.86, RPD = 2.68) were statistically equivalent [50]. These results suggest that in situ proximal sensing under varying soil moisture conditions combined with ML can achieve similar good soil predictions as those derived in controlled conditions in the laboratory.

ML algorithms have been widely applied in the emerging field of predictive soil modeling using portable spectroradiometers that characterize soils under field conditions. Portable VNIR and MIR approaches have shown promising results using PLSR modeling of soil carbon and other soil health properties when compared to lab-based diffuse reflectance spectral measurements [27, 51,52,53,54,55,56,57]. According to Hutengs et al. [55], portable VNIR and MIR instruments provided accurate models of various soil physicochemical properties (R2 between 0.72 to 0.99) that showed some influence by the soil moisture state (dry vs. field-moist). Validation models for SOC achieved an R2 of 0.82 (dried, VNIR), 0.88 (dried, MIR), 0.57 (field-moist, VNIR), and 0.72 (field-moist, MIR). In the study presented by Semella et al. [56], SOC predictions from both VNIR and MIR spectra collected with portable spectroradiometers were equally highly reproducible on average with slightly higher robustness in MIR. Results showed that the contributions of spectral variation (∆RMSE < 0.4 g kg−1; RMSE: Root Mean Square Error) and reference SOC uncertainty (∆RMSE < 0.3 g kg−1) to spectral modeling errors were small compared to the difference between the VNIR and MIR spectral ranges (∆RMSE ~1.4 g kg−1 in favor of MIR). Studies with handheld single-sensor instruments, such as ASD Labspec 2500 [51], Quick Carbon Reflectometer [58], Agilent 4300 handheld FTIR [53], AgriSpec [57] (Sharififar et al., 2019), NeoSpectra [27,57,59], NanoQuest [60], and Hamamatsu C12880MA [27] demonstrate the capabilities to sense SOC and other soil properties, though with variable results based on the sensors capabilities. One major disadvantage is that these portable instruments require soil samples to be extracted to be sensed in the field and they do not allow in situ continuous sensing along soil profiles. These kinds of quasi-in situ VNIR sensing systems require soil cores to be first extracted and then scanned using a field spectroradiometer [44]. Tractor or truck-mounted sensors cover the full VNIR spectral range, but due to the vehicle movement during data collection, often, the uncertainty in soil predictions can be substantial [27,61]. Soil sensors that do not possess in situ penetration capabilities severely limit characterization of soil spatial variability, especially in crops with extensive rooting system.

A comprehensive characterization of a suite of soil health and other profile properties and attributes in agriculture applicable to a wide variety of cropping systems (e.g., specialty crops, row crops, different crop species) calls for multiple sensors to be used in combination that are fully integrated into a soil sensing system. Often single sensor instruments are applied separately to map specific soil characteristics and then data are fused later during the data processing and modeling phase [62]. For example, individual sensors such as apparent electrical conductivity (ECa) to map soil salinity [63], portable X-ray fluorescence (pXRF) spectrometry for elemental and soil fertility characterization [64], and high capacity tensiometers, microwave-based approaches, and others for soil moisture sensing provide specialized applications. Schmidinger et al. [65] compared the model performance of six independent in situ proximal soil sensors, one remote sensor (Sentinel-2), and all sensor data fused together to predict SOM, phosphorus (P), magnesium (Mg), potassium (K), moisture, and pH with multiple ML algorithms. Five out of six soil properties achieved an R2 ≥ 0.80 often with various combinations of individual sensors, while unsurprisingly the improvement derived from fusing an increasing number of sensors was subject to diminishing return. Similar testing of soil model performances to assess the effectivity of multiple single-sensor combinations (less than max. of 4) and fused sensor data were presented by Chen et al. [66] (2021), Tavares et al. [67], and Xu et al. [68]. Vasques et al. [69] applied multiple sensors (ECa, apparent magnetic susceptibility meter, gamma-ray spectrometer, water content reflectometer, cone penetrometer, and pXRF) in a pasture field and found that multiple soil sensor data fused together improved soil predictions for all soil properties relative to single sensors. The pXRF data produced the best predictions for SOC, clay content, and BD, 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, CEC, and soil volumetric moisture, based on independent validation. These findings suggest that different combinations of sensors are needed to provide inference on a variety of soil physical and chemical properties.

Although the integration of multiple sensors into a mobile platform has sparked profound interest in the agronomic and soil science communities, fully integrated systems are rare and typically limited to few sensors. An early attempt at a multi-sensor system for soil physical properties was presented by Yurui et al. [70]. The Veris P4000 multi-sensor instrument can collect VNIR spectra, ECa, and cone index (CI) penetrometer readings up to 1 m depth. In Pei’s study in two fields in central Missouri, U.S. the Veris P4000 achieved modest results in cross-validation mode with average R2 values across all soil properties (SOC; total nitrogen, TN; soil texture, CEC, Ca, Mg, K, and pH) for PLSR, neural network (NN), Regression Trees (RT), and RF were 0.59, 0.46, 0.39, and 0.45, respectively. Though few properties achieved promising results with PLSR (e.g., R2 of 0.81 for SOC), while some properties showed weak model fit (R2 of 0.37 for sand content). A multi-sensor robotic platform with a modular sensing box that includes VNIR, thermal camera, two visual cameras forming a stereo couple, and an Inertial Measurement Unit (IMU) that provides navigational data mounted on an autonomous vehicle to generate 3D ground maps for precision agriculture applications was described by Milella et al. [71]. Other multi-sensor soil systems are static and intended for real-time sensing at only one specific location. For example, a buried soil probe containing electrochemical sensors in a hygroscopic membrane to monitor soil nutrient concentrations in real-time was combined with an air probe that collects information regarding environmental conditions, gaseous emissions (esp. NH3, N2O, CH4) just above the ground, and smart data loggers connecting to the Internet of Things (IoT) cloud [72] (Balan et al., 2020). Such static soil sensor systems lack mobility to collect data across farms and cropping systems to optimize climate-smart and practical agricultural management.

In this paper we present research using an in situ proximal soil sensing system designed and deployed by LandScan, LLC (Davis, CA) that includes a multi-sensor probe, software, and equipment to deploy (DSC System). The research objectives include investigating: the capabilities of the DSC System to predict various soil health and management-related properties as well as directly predicting crop metrics without the use of ex situ soil samples and laboratory analytics; and the effects of mismatched scale between high-resolution in situ proximal sensor data and coarser resolution ex situ soil laboratory measurements to develop soil and plant prediction models used to create a digital twin. We critically discuss the limitations of the contemporary paradigm to ground-truth soil sensor data with laboratory-based ex situ soil measurements and present an alternative method that focuses on measured soil-crop responses.

Study Area

Data collection for this study was conducted across three almond management blocks located on commercial ranches in central California (Figure 2). The first ranch is positioned near the San Joaquin River, southwest of Madera in Madera County, while the remaining two are in Kern County, southwest of Bakersfield. Detailed descriptions of the location, size, crops, soils, climate, are found in Table 1. The almond trees were 7 to 12 years old and irrigated using drip or micro-sprinkler irrigation. The Central Valley of California has a Mediterranean climate, characterized by hot, dry summers and cool, wet winters. Trees are planted on linear berms that extend 10 to 20 cm above the lanes. The lanes have a cover crop in the winter/spring but are typically cut back in mid-summer to facilitate ground preparation for harvest in late summer. The berms are kept free of cover crops on these sites.

The block KG 18-19 (size: 35 ha) is north of the San Joaquin River, SSR 35-1 (size: 25 ha) and ST-15 (size: 31 ha) are located adjacent to canals. ST-15 previously had a drainage or canal running through it and was previously part of the adjacent cattle ranch. ST-15 is split into two parts by a gravel ranch road.

2. Materials and Methods

2.1. Digital Soil Core System and Probe

Our research employed the DSC System that includes the integrated components of the DSC Probe (Figure 3), software, and equipment to deploy. The DSC Probe is a multi-sensor probe that includes 1) tip stress, 2) sleeve friction, 3) dielectric permittivity, 4) electrical conductivity, 5) Micro-Electro-Mechanical System (MEMS) microphone, 6) video microscope, and 7) visible and near-infrared (VNIR) diffuse reflectance spectrometers [73]. The DSC probe can penetrate the soil up to 120 cm in this configuration. Tip and sleeve stress measurements are indicators of soil strength [74] which is spatially and temporally variable. The DSC Probe incorporates a 60-degree, 1-inch diameter conical tip [5] (Rooney and Lowery, 2000). A pair of steel electrodes in the tip of the DSC Probe are separated by an insulating element and used to measure volumetric water content and electrical conductivity via rapidly multiplexed measurements of direct current (DC) electrical resistance and apparent dielectric permittivity at a frequency above 50 MHz from which water content is inferred. The dielectric permittivity of the soil is recovered via calibration to known standards and converted to volumetric water content (VWC) using well-established relationships [75,76,77]. An embedded microelectromechanical systems (MEMS) digital microphone records acoustic emissions produced by the penetration process, as soil particles are scraped and re-arranged due to penetration displacement [78,79]. The sound is affected by soil texture and structure, compaction state, and water content, making the microphone sensitive to several important soil parameters. Two sapphire windows permit video microscope imagery and VNIR DRS with optics and lighting optimized for subsurface microscopy at penetration speed. Uniform, consistent illumination is synchronized to the video frame rate. The videos were captured using the Advanced Video Coding (AVC), H.264, video compression standard. The video is captured in H.264, Red Green Blue (RGB) frames and extracted for processing. The microscope produces RGB color imagery (2.3 x 1.2 mm) with 1-μm pixel resolution and a spatial density of about 15 images per cm with 50% overlap of adjacent images. Optical resolution of 3 μm is confirmed using a MIL-STD-150A resolution calculator (#38-257, Edmund Optics). VNIR DRS data is conducted at a rate of four scans per second, with a push rate of 2 cm/second resulting in approximately 2 VNIR readings per cm. The downhole optical design and proprietary optical fiber bundle of the VNIR system is optimized for maximum signal-to-noise ratio (SNR) in spectra collected by spectrometers located above ground and external to the DSC Probe from Ocean Optics (QEPro and NIRQuest) in a custom enclosure engineered for environmental protection and precise thermal control. The QEPro has a spectral range of 350 to 950 nm and full-width half-max (FWHM) optical resolution of 1.2 - 6.87 nm. The NIRQuest has a spectral range of 900 to 2,500 nm and full-width half-max (FWHM) optical resolution of 6.3 nm. The DSC System includes a string potentiometer used to register the DSC Probe depth during penetration

In contrast to conventional core retrieval and laboratory analysis, the DSC Probe data collection method preserves vertical spatial variability, differentiates thin layering, and accurately references soil parameters to depth. Other advantages that in situ proximal sensing can provide over traditional ex situ soil coring, compositing, and homogenization of soil samples include observing the in situ distribution of soil water within the structural arrangement revealed, and the in situ bulk electrical conductivity rather than that of saturated paste extract for example.

The integration of multiple independent proximal soil sensors in the DSC Probe enhances the capacity to capture a comprehensive picture of soil properties and the in situ relationship to each other. Each sensor modality offers a unique perspective on soil properties, and when combined, they provide a multifaceted characterization of the soil profile (example sensor vertical plots are provided in Figure 4 and imagery in Figure 5). An important consideration in the development of the DSC System is the interplay of orthogonality of sensor modalities and degrees of freedom in the sensor data and soil parameters of interest. Accounting for the dimensionality of information within individual sensors such as video, audio, and spectrometry, the DSC System provides over 1200 sensor output values for each cm of soil it encounters.

2.2.1. Soil Data Collection

The in situ and ex situ soil data collection took place in October and December 2023. Figure 6 shows the DSC System in operation in an almond orchard.

DSC sampling locations and collocated soil cores were targeted using c-means clustering [79] applied to the EM data. The c-means clustering algorithm was used to find 6 clusters and identify one DSC target location per cluster. Additional DSC observations were obtained for commercial mapping purposes but were not included in this study.

Both DSC sensor measurement profiles and physical soil cores were obtained in triplicate at each target location. All were acquired from within an area measuring approximately 1-m by 1-m at each target location (Figure 7), between the center and the shoulder of the berm on the tree-row berm between two almond trees.

Each DSC Probe measurement profile extended to about 1.2 m below ground. To assure the highest accuracy of spectral reflectance data, free of instrument thermal drift and other factors, the DSC System automatically performs a series of reference dark current scans at the terminus of every digital profile. With the sapphire window embedded more than 1 m deep in the ground, free from any possibility of ambient light, the illumination source is shuttered, and a dark current reference measurement is obtained. Conversion of raw spectral scan data to reflectance spectra considers the nearest-in-time dark current reference scan along with the nearest-in-time white reference scan obtained by covering the sapphire window with a Spectralon® diffuse reflectance standard (Edmund Optics Stock #54-302) and triggering the control software to acquire a series of reference reflectance scans. Processing of DSC sensor data is described later.

Electromagnetic induction (EM) data was collected along rows in the almond orchards to help understand soil variability patterns with a Dualem-1HS (2 distances x 2 orientations) giving 4 channels of apparent electromagnetic conductivity to 4 depths of exploration (30, 50, 80, 160 cm). The EM was driven down each row of the mapped orchards. A Real-Time Kinematic (RTK) GPS was used for georeferencing of EM data. The data were then processed on-the-fly using LandScan data collection software to remove physical and temporal offsets between GPS and EM and vertical offset between GPS and the ground. Results were filtered using a windowed standard deviation filter and interpolated to rasters using a thin plate spline algorithm (Minimum Curvature interpolation algorithm in Datum Workstation a geospatial analysis system formally known as TNTmips, LandScan, 2023).

Physical cores were obtained using a 122-cm (48-in) core barrel, with plastic liners having an inner diameter of 41 mm (1.6 in). The soil cores from each location were aligned in a tray, starting at the top and extracted from the plastic tube. Any obvious horizon breaks in the soil were aligned between the cores. Cores were then broken into 10-cm (6-in) horizons across all three cores. Soil that appeared consistent for each horizon (at least volumetric equivalent to ~500g) was bagged for lab analysis. 1 bag per horizon across 3 cores. Three cores were used to keep horizon thickness small while providing sufficient sample volume to the lab. Samples beyond 110 cm (42 in) were not sent to the lab. Soil samples were labeled with DSC push identification numbers (IDs) so that they could be matched to the DSC sensor data for training.

In total 60 soil cores and 60 DSC digital profiles were collected within the whole study area. Refer to Table 1 for the number of soil samples submitted by study area. Approximately 6-8 sampling depths were selected from each ex situ soil core location and sent to a commercial laboratory (Dellavalle Laboratory Inc, Fresno California) for analysis. Soil analytical measurements included organic matter (OM, loss on ignition), particle size (sand, silt, clay measured by the hydrometer method) and a complete soil fertility package. Out of the measurements in the soil fertility package, Boron (B), Calcium (Ca), Copper (Cu), Zinc (Zn), and pH were evaluated in this study. References to the laboratory methods used are included in Table 2. Note that Nitrogen (N) was excluded because the concentrations at each site were very low with no data distribution to measure against.

2.2.2. Crop Data Collection

For this study, the almond crop vegetation was characterized utilizing the Digital Vegetation Signature ^TM (DVS) technology developed by LandScan [80]. Each site was flown mid-season for the study with a DGI Mavic M3M multispectral unmanned arial vehicle (UAV) at an altitude of 120 m. The UAV has an RGB camera, a multispectral camera and a built in GPS. The imagery was processed using the Rig Camera Alignment tool in Datum Workstation. A spectral calibration was performed against ground control targets prior to mosaicking in Datum Workstation which was then used to produce the final orthorectified mosaics. The mosaics were processed into a vegetation vigor index (VVI), a pigmentation index (PI), and numerous other indices using proprietary algorithms in Datum Workstation. The richness of the combination of both spectral and spatial data combined reveals many new features in and about the data that provide valuable input to future analytical processes and integration into the LandScan Digital Twin for Agriculture [81].

In addition to the ortho rectified imagery, the orthorectification process also results in a digital surface model (DSM). The digital surface model was used in conjunction with a digital terrain model (DTM) acquired from the U.S. Geological Survey (USGS) National Elevation Dataset program along with a vegetation raster to create a vegetation height raster. The vegetation height raster is used to determine the location of each tree in each orchard block, give them an identifier and establish various canopy masks. One example is that for each pixel in the canopy mask, height is multiplied by VVI and then summed to form a total Crop Productivity Index (CPI). This approximates canopy volume and density, or total canopy biomass (e.g., volume x density should equal mass), which relates to the fractional amount of photosynthetically active radiation (fPAR) that can be absorbed by each tree. In theory, an almond tree’s productive capacity is limited by fPAR [82]. The DVS data collection process results in a finite number of ‘named’ indices, but also produces many new data features and relationships that empower deeper learning opportunities for advanced analytics (Figure 8). Many of these features and relationships are integrated into the digital twin for exploratory and discovery purposes as additional agronomic metrics become available. These data enable a wide range of opportunities to advance and improve on the approach taken in Method D in this study.

2.3. Data Pre-Processing and Harmonization

The multi-sensor data that were collected with the DSC System were screened to identify outliers, noise, and missing data as part of quality assurance procedures in LandScan’s DSC data collection software. Data collected from all three profiles were then merged spatially to obtain one representative digital profile to integrate with soil properties and crop responses for modeling purposes.

2.4. Spectral Data Processing

Using reference standards and dark current measurements as described above, spectral reflectance from the DSC Probe was first computed. In this work, reflectance spectra were converted to absorbance spectra before applying Standard Normal Variate (SNV) Transformation and Savitzky-Golay (SG) filtering [83]. The SNV minimizes multiplicative effects such as baseline shifts and light scattering in spectroscopic data [84]. SG filtering was applied to remove noise and improve the signal-to-noise ratio of spectral data while preserving spectral features. For this purpose, 1st differential order and 2nd polynomial order with 11 window sizes were used. Spectral pre-processing was performed using the prospectr package in R (https:// CRAN.R-project.org/package=prospectr).

2.5. Processing of Digital Soil Images

Image color metrics mean hue, value, and saturation (HSV) were extracted from each DSC Probe microscope image, as well succolarity, a metric of image structure (de Melo et al. 2008), for consideration in the analysis. Of these image metrics, succolarity curve difference, color saturation, color hue and color value were found to have significance in the final ML model. Succolarity was originally developed to measure the flow of water through canal systems from satellite images [85] with additional flow-related applications suggested by de Melo and Conci [86]. LandScan applies succolarity algorithms to quantify the potential for percolation flow through porous media in an image in the analytics software. The determination of succolarity begins with binary masking of the image based on a threshold value below which a pixel is considered to represent void (pore space), and above which a pixel is considered to represent structure (soil matrix). The binary image is then flooded with a theoretical ‘fluid’ from each of the four edges of the image boundary, and the proportion of the total image penetrated by the fluid from each direction of flooding is computed. The four values are then averaged into a single succolarity value. This approach to computing succolarity, by Leavitt et al [87](2021), approximates the methods explained in de Melo and de Melo and Conci [85,86]. Since the succolarity value thus computed is a function of the threshold chosen for the binary masking operation, we generated multiple values of succolarity as a function of the threshold value chosen which comprise a succolarity curve. This curve tends to exhibit a sigmoidal shape and the metric we call succolarity curve difference is the normalized difference in image masking thresholds between the start and end of the rise in the succolarity curve [87].

2.6. Processing of Audio Data

Audio data from the MEMs microphone was recorded in Waveform Audio Format (WAV) for storage. The WAV file was processed in Python by converting to a numpy array and running a 3-kHz high pass Fourier transform filter, followed by binning into five bins of 4-kHz bandwidth, a sound pressure level for band and the total sound pressure level.

2.7. Processing of Other Sensor Data

All DSC System sensor-derived data were harmonized to co-registered 1-cm depth increments in LandScan DSC processing software. The outputs of all DSC Probe sensors and the DSC System string potentiometer were used to register the DSC Probe depth during penetration and are each associated with a time stamp during data acquisition. Because the sensors and their contact with the soil each occupy a different position along the DSC Probe as it advances through the vertical profile, each increment of soil is encountered by a different sensor at a slightly different time. To co-register the readings from all sensors with respect to depth, given that slight variations in penetration speed may occur during acquisition of a sensor profile, each time series of sensor readings is first independently indexed to depth and then re-sampled relative to a common index of equally spaced depth intervals, such as every 1-cm. Depth co-registration is achieved by applying a sensor-specific depth offset to each sensor in the probe based on its relative position in the DSC Probe, then computing the depth each sensor was at when each of its readings was recorded, then re-sampling the readings from each sensor independently using cubic spline interpolation to conform to a uniformly spaced set of depth values distributed over the depth of the profile with a depth referenced to zero depth at the ground surface.

2.8. Data Feature Selection

The Boruta feature selection algorithm was applied to reduce the dimensionality of the massive data cube of sensor data by identifying the most relevant sensor output for predicting soil properties (Method A, B and C) and crop responses (Method D). It is one of the widely used variable selection methods in soil spectroscopy to deal with the multi-collinearity of data [88,89,90]. Boruta trains an RF model using a combined dataset of original and shuffled features and evaluates the variable importance (Z score) for each predictor. Then it checks whether a real predictor has higher importance (RMSE) than the best of its shadow predictors to decide important and unimportant features. In this study, all the high-resolution (1 cm) sensor data were used as features data in the RF classifier from the Scikit-Learn library in Python to select the important features for individual soil properties and crop responses [91].

2.9. Comparison of Training Methods

Four different methods were used to assess the model performance of soil health and management properties. Modeling was performed with PLSR using the leave-one-out cross-validation [92]. The goal was to determine the best method of assessment between A, B, and C as compared to the laboratory, then use that method to compare to Method D in predicting crop response (Figure 9).

Method A reduced the high-resolution DSC Probe sensor data to discrete values paired with layer-based soil laboratory measurements. All the high-resolution (1-cm) sensor data was averaged to match the length of segments of ex situ soil cores sent for the laboratory analysis. In essence, for each laboratory measurement, one array of DSC Probe sensor data was used in the calibration models. We used the leave-one-out cross-validation method and PLSR on the sensor and soil analytical data from each of the 15-cm layer increments for all cores.

Method B used stochastic distributions of the DSC Probe sensor data paired with layer-based ex situ soil laboratory measurements. Here stochastic distributions of all DSC Probe sensor data in the model were used for PLSR modelling. In this method, soil analytical laboratory data were matched with the minimum, maximum, standard deviation, and mean sensor data associated within a 15-cm layer. For validation, the arrays of minimum, maximum, standard deviation, and mean DSC Probe sensor data for corresponding soil analytical laboratory samples were used in leave-one-out cross validation.

Method C allocated the same soil analytical laboratory data to each one of the high-resolution multi-sensor data within a layer. PLSR models used all high-resolution (1-cm) DSC System data corresponding to soil analytical laboratory measurements. Since laboratory measurements were only available for each layer, the same laboratory data values were matched with each corresponding high-resolution DSC System data. For validation, however, we ensured that for each laboratory sample left out during cross validation, every high-resolution DSC System data increment corresponding to the sample that was left out was also left out. The predicted soil properties were averaged for each 15-cm layer.

Method D linked the high-density in situ DSC System data directly to DVS crop responses (crop performance and behavior metrics), bypassing costly laboratory soil analysis. In this approach, DVS crop responses such as the Crop Productivity Index (CPI), Canopy Area, and Canopy Volume were directly predicted from the DSC System to avoid the laboratory measurements of soil properties. Since crop responses are single measurements of each location, the optimum depth for aggregating sensor data was determined. A few different soil depth intervals (0-20, 0-30, 0-60-cm) were considered to find out the optimum depth of the DSC System data that predicts crop response with comparatively higher accuracies. Finally, based on soil health and nutrient management opportunities in almonds the 0-30 cm depth was considered for the analysis. The same crop response was matched with the array of high-resolution (1-cm) DSC System data for training and validating purposes. Then, the predicted crop responses were averaged and compared with the observed DVS CPI, Canopy Areas, and Canopy Volumes.

2.10. Modeling Approach

The PLSR with leave-one-out cross-validation approach was used to estimate soil properties (Method-A, B, and C) and crop responses (Method-D). All the samples from individual fields were used for training the calibration model except one sample that was used for validating the calibration model. The number of components (n-component) used to obtain the lowest RMSE between the measured and estimated soil properties and crop responses in the training model, was used for validation purposes. Modeling was performed using the Python programming language with the “PLSRegression” from scikit-learn 1.2.1 package. For evaluating the performances of all four Methods, R2, RMSE, RPIQ, and bias of the modeling were used.

3. Results

This section may be divided by subheadings. It should provide a concise and precise description of the experimental results, their interpretation, as well as the experimental conclusions that can be drawn.

3.1. Feature Selections for Modeling

Before training the predictive models for the estimation of soil properties, the Boruta feature selection algorithm was applied to over 1,200 features per cm in the DSC System dataset to obtain the importance of sensor data that will be used in the models. The top 20 important features of soil properties and crop responses are shown in Figure 10 for Zinc (Zn) and Figure 11 for CPI. The other Boruta graphs can be found in Appendix A. For both the soil properties and crop responses, all the features that had a Z score of more than the maximum shadow value were used for the predictions. The range of features used in soil property predictions in this study varied from 36 to 83.

For almost all soil properties, DSC Probe sleeve friction, penetration resistance, friction ratio and 1600-2000 nm wavelengths from VNIR spectra were found to be the most important features. From the whole VNIR spectral range, near-infrared (NIR) wavelengths seem to be important for all the soil properties except for Ca.

For crop responses, color (hue), color (saturation), color (value), soil moisture, penetration resistance, succolarity curve difference (SCD), electrical resistance, macro porosity, sleeve friction, and friction ratio were the important features along with several bands from VNIR spectra. Similar to feature importance for soil properties, the NIR region has more important features compared to the visible region of the spectra for crop response.

3.2. Predictive Accuracy of Soil Properties Modeling Methods

This study compared prediction capabilities of three different methods for modeling various soil properties (OM, Sand, Clay, Silt, B, Ca, Cu, Zn, and pH) using DSC digital soil profiles. The results for R2, RMSE, Bias, and RPIQ are reported in Table 3.

Method C had the highest R2, the lowest RMSE, and highest RPIQ across all 9 soil properties at all 3 sites. It also had the highest mean R2, lowest mean RMSE, and highest mean RPIQ across the 3 sites for all 9 soil properties, indicating that Method C leads to more accurate and robust results, than the other two methods. Method B only slightly outperformed Method A for some properties (OM, Sand, Clay, Silt, pH) but not all properties (B, Ca, Cu, Zn).

3.3. In situ DSC System to Ex situ Laboratory Properties to Digital Crop Performance vs. DSC System to DVS Digital Crop Performance.

We compared the performance of two models that predict crop performance from in situ and ex situ soil data. In the first approach, we used the best performing method of predicting soil properties, Method C, and predicted the soil properties, and then used the soil properties to predict CPI, canopy area, and canopy volume as measured by a UAV across all three sites.

In the second approach we directly predicted CPI, canopy area, and canopy volume based on the in situ DSC System data directly, without predicting ex situ laboratory soil property values (Method D). For CPI, Method D had an R2 between 0.72 to 0.75 and an RPIQ between 1.13 to 1.64 (Table 4), whereas the prediction of CPI using Method C had an R2 between 0.54 to 0.67 and RPIQ between 0.63 and 0.85. Method D had a higher R2 and lower RMSE than Method C, indicating that in situ DSC System data (input data) to digital crop response models show higher accuracy than a more complex approach that sequentially models in situ DSC System data (input data) → soil properties → crop response models.

4. Discussion

The results show that the sensor data collected from the in situ DSC System has the potential for estimating soil properties and crop responses with the support of chemometrics modeling. All three methods (Method A, B & C) used in this study showed some correlation between the DSC System data and various soil properties, but Method C exhibits the highest prediction accuracies compared to other methods. The RPIQ for Method C ranged between 1.66 (Zn) to 2.94 (pH) which are compelling results for soil predictions using in situ DSC System data. The R2 for soil properties (Method C) ranged between 0.61 (Cu) up to 0.79 (pH); and 0.73 (OM) suggesting improved results compared to other proximal soil sensor applications. For example, models to predict SOC stock using field VNIR spectral data in a study in France achieved an R2 between 0.52 to 0.86 and RPIQ between 1.61 to 4.49 in validation mode [51]. On experimental plots in the Canadian provinces with humid soil moisture regime SOC concentrations modeled by VNIR spectra achieved an R2 or 0.54 (MIR) and R2 of 0.49 (VNIR) [52]. In a study in Germany, topsoil SOC (%) predictions (validation mode) using Veris full-range VNIR device and PLSR modeling showed modest R2 (0.55) and RPIQ (2.05), the Hamamatsu sensor showed poor performance with R2 of 0.29 and RPIQ of 1.67, while the NeoSpectra results were slightly better with R2 of 0.48 and RPIQ of 2.00 [27]. Many soil sensor applications focus only on predictions of SOC or OM, while the DSC System presented in this study has a much broader range to model a suite of different soil properties.

One of the major advantages of the DSC System is the high resolution (< 1 cm) of data acquisition from all the sensors. In Method A, the mean of all 1-cm data in a horizon was used to train and validate the model. For Method B, instead of only using mean data, minimum, maximum, standard deviation, and mean of all sensor data in a soil layer were used. However, in Method C, the full potential of the high-resolution DSC System data was used, which allowed training of the model retaining the variation of sensor data along soil profiles. Higher prediction accuracies from Method C for all soil properties indicate the importance of recording high-resolution sensor data for accurate soil predictions.

From Boruta feature selection (Figure 8), it was observed that most of the important features were obtained from VNIR spectral ranges specifically for soil textures (Clay, Silt and Sand), OM and pH. Prediction accuracies for these properties are better compared to the other soil properties estimated using the DSC System. Other DSC Probe data, such as sleeve friction, penetration force, friction ratio and color saturation, are the most important features, as these appeared in most of the Boruta important feature plots. Several features, such as succolarity curve difference, electrical resistance, macro porosity and soil moisture, were also found to be important for the estimation of pH and Zn. Most of the important features were found from the VNIR spectral ranges for nutrients except Zn. For the crop responses all the color properties (hue, saturation, and value), penetration resistance, sleeve friction, friction ratio, soil moisture and micro porosity were found to be the most important features followed by VNIR spectra.

It is practically not possible to obtain an objective, spatially accurate map if ex situ soil sampling is utilized on large farming operations (large fields and ranches). Conventional wet chemistry analyses involve extracting soil cores from the field, transporting samples and processing samples for laboratory analyses. Maintaining all these standard protocols disturbs the original condition of soil samples [93,94]. This study brings into serious question the accuracy and applicability of conventional ex situ soil sampling and laboratory practices for advanced agronomic analytics and negates the opportunity to produce a digital twin. Many human (e.g., handling of soil samples or cores) and laboratory measurement errors may occur without even acknowledging and quantifying them explicitly.

This study shows that multi-sensor data collected using the DSC System can rapidly and objectively estimate multiple soil properties. All these data were collected in situ and within a fraction of the time for extracting ex situ soil cores in the fields and tested in the laboratory using wet chemistry analyses. Importantly, the DSC System reduces the time and cost of characterizing the soil profile by reducing or even omitting the expenses for extracting ex situ soil cores, processing, transporting to commercial soil testing laboratories and wet chemistry analyses.

Since the density of the DSC System data (< 1 cm scale) is substantially higher than discrete sample extraction in different soil layers by traditional soil analysis in the laboratory one may argue that the real “gold” standard are the DSC System data. These sensor data are collected in close proximity to the soil matrix under actual field conditions providing a more direct way to characterize soils than conventional ex situ soil surveys. Therefore, our study lays the foundation to shift the paradigm of future soil sensor applications to focus directly on sensor data (e.g., VNIR hyperspectra, porosity derived from digital micro-images) and crop responses (e.g., canopy density) rather than soil interpretations (e.g., OM or soil texture).

In this study, DSC System data was successfully used to estimate crop responses in all three fields (Method D). Bypassing the estimation of the soil properties to estimate the crop responses directly from the sensor data can potentially offer a more streamlined, objective, efficient, and accurate approach to precision agriculture. This approach avoids the risk of errors associated with indirect, subjective, and analog measurements of some soil properties that are not related to the sensor data used for modeling. Directly estimating crop responses from soil sensor data simplifies the data processing pipeline by eliminating intermediate steps (e.g., soil property predictions that match laboratory measurement methods), leading to faster and more efficient data analysis that can be used for precision management of crops.

5. Conclusions

We demonstrated that the DSC System makes robust soil property predictions across multiple soil properties and study areas in central California using a standard machine learning approach. Use of the technology in other crops and soil regions and applications of advanced ML algorithms will further improve sensor-driven soil and crop modeling approaches that promise substantial future cost-savings. The innovative DSC System facilitates collection of standardized soil signatures from multiple concurrent sensing modalities spatially co-registered within specific soil profiles.

The integration of the DSC System’s multiple sensing modalities better conditions the ML model solutions that infer specific soil properties or soil-crop relationships from the sensor data. This provides predictive performance that is superior to sensor systems with fewer modes such as proximal soil sensors that use diffuse reflectance spectroscopy alone or with fewer complements. Each additional sensor modality added confidence. In addition, integrated multi-sensor data collection from a single device is less error prone than multiple single sensor systems used in combination, which suffer from sensor displacement and disharmonious sensor resolutions that require more extensive post-processing of data to correct and more extensive soil sampling to support.

The best method of training a model to predict soil laboratory data was Method C, which involves treating each centimeter of sensor data as a separate measurement that is paired with a soil lab measurement. Our results suggest that the machine learning algorithms can learn more from a high density of sensor data retaining the spatial variation in soil characteristics along a soil profile compared to aggregating sensor data to coarser scale (i.e., collecting a soil sample representing a soil layer that is then analyzed in the lab) to match conventional soil surveys. The robust model performance underpins the importance of vertical scale when characterizing soil properties with multiple sensors outperforming traditional soil surveys. The collection of in situ sensor data in soils is a prerequisite to create realistic digital soil twins that cannot be achieved with soil core extraction and/or through conventional soil laboratory analysis. Thus, we envision a new technology informed “gold” standard for digital soil mapping employing a multi-sensor in situ proximal sensor suite combined with AI modeling rather than the traditional standard of discrete soil sample extraction and ex situ analysis.

The best method of training a model to predict crop productivity was Method D (sensor data → crop responses) which outperformed the more complex approach using sequential modeling (sensor data → soil properties → crop responses). These results suggest that direct sensor-crop modeling has less error and higher accuracies than sensor-soil-crop modeling which suffers from error propagation lowering overall model performance. From a statistical perspective, clearly the path of sensor data → crop responses modeling is preferable. The analytical potential of combining a full DSC digital soil profile as a source of calibration for drones, airplanes, and satellite data is compelling.

While direct prediction of crop productivity could be useful in determining productivity potential and inform and improve certain agronomic practices, soil property prediction will still be valuable for describing the below-ground factors that affect that potential. In turn this helps determine what actions a grower can take to improve crop productivity in their fields such as the production of variable rate (VR) fertility and soil amendment maps or adapting irrigation practices to optimize based on soil variability

Multi-modal in situ proximal soil sensing systems such as the DSC System present immense potential to transform soil-crop digital mapping and modeling. We continue to acquire DSC System data and corresponding ex situ soil cores for laboratory analysis from numerous locations in California and other locations (e.g., Australia). All the samples with ground truth data are used to develop a DSC data library that is currently used to train ML models using data that comprises large variations in soil properties and conditions (e.g., moisture content, fertility, etc.). Over time, we plan to expand the DSC data library with ex situ soil samples and measurements across the U.S. and major agricultural areas world-wide. The goal is to reduce and eventually eliminate the need to extract soil cores from every field. A sufficient quantity of data will enable the application of more data-hungry, deep learning models that will use the diverse and extensive dataset for training the prediction models more efficiently as new soils are added to the library. The DSC data library will enable exploitation of the full potential of the DSC System’s speed, cost, and reproducibility advantages by estimating soil properties and crop responses from any field in the future using only DSC data. These developments are essential to inform decision support systems that truly optimize climate-smart agricultural management, high-accuracy soil carbon accounting, precision agriculture applications, and the installation of management unit level digital twins. Looking forward, this high spatial and information density data cube will be the type of input necessary to run quantum computing models for future agricultural decision support particularly in intensively managed cropping systems facing resource constraints.