Total nitrogen estimation in agricultural soils via aerial multispectral imaging and LIBS

Measuring soil health indicators (SHIs), particularly soil total nitrogen (TN), is an important and challenging task that affects farmers’ decisions on timing, placement, and quantity of fertilizers applied in the farms. Most existing methods to measure SHIs are in-lab wet chemistry or spectroscopy-based methods, which require significant human input and effort, time-consuming, costly, and are low-throughput in nature. To address this challenge, we develop an artificial intelligence (AI)-driven near real-time unmanned aerial vehicle (UAV)-based multispectral sensing solution (UMS) to estimate soil TN in an agricultural farm. TN is an important macro-nutrient or SHI that directly affects the crop health. Accurate prediction of soil TN can significantly increase crop yield through informed decision making on the timing of seed planting, and fertilizer quantity and timing. The ground-truth data required to train the AI approaches is generated via laser-induced breakdown spectroscopy (LIBS), which can be readily used to characterize soil samples, providing rapid chemical analysis of the samples and their constituents (e.g., nitrogen, potassium, phosphorus, calcium). Although LIBS was previously applied for soil nutrient detection, there is no existing study on the integration of LIBS with UAV multispectral imaging and AI. We train two machine learning (ML) models including multi-layer perceptron regression and support vector regression to predict the soil nitrogen using a suite of data classes including multispectral characteristics of the soil and crops in red (R), near-infrared, and green (G) spectral bands, computed vegetation indices (NDVI), and environmental variables including air temperature and relative humidity (RH). To generate the ground-truth data or the training data for the machine learning models, we determine the N spectrum of the soil samples (collected from a farm) using LIBS and develop a calibration model using the correlation between actual TN of the soil samples and the maximum intensity of N spectrum. In addition, we extract the features from the multispectral images captured while the UAV follows an autonomous flight plan, at different growth stages of the crops. The ML model’s performance is tested on a fixed configuration space for the hyper-parameters using various hyper-parameter optimization techniques at three different wavelengths of the N spectrum.

Soil health indicators are a composite set of measurable physical, chemical and biological properties which can be used to determine soil health status. Among the chemical indicators, we particularly focus on nitrogen (N) because N is the most limiting nutrient in many of the world's agricultural areas 1 . Insufficient use of N causes economic loss, in contrast, excessive use of N implies wasting fertilizer, causes nitrate pollution, and increases the cost 2,3 . Nitrogen treatment can account for up to 30% of the total production cost 4 .
Chlorophyll meter (CM) measures the chlorophyll content of crops to estimate their N nutrition status. In recent years, the use of CM has increased among researchers and farmers 5,6 . For instance, N application rates for corn were determined using the adjusted R 2 of the relationship between nitrogen rate difference (ND) and CM readings 7 . However, CM-based methods fail to capture the spatial variability that is often present within the field. For N management, determination of spatial patterns is necessary but requires collection and analysis of a large number of samples which is labor-intensive and time-consuming 2,6 .
Satellite-based remote sensing is one alternative to ground-based measurements. Satellite-based techniques utilize images at the spectral level for crop growth monitoring and real-time management [8][9][10] . For instance, vegetation indices (VIs), evaluated using the data obtained from satellite-based multispectral sensors, have been used to detect the N stress at V4-V7 (4-7 leaves with visible leaf collar) stages [11][12][13] . However, satellite-based sensing www.nature.com/scientificreports/ suffers from lower spatial and temporal resolution, and sensing disruption may occur during image acquisition in some areas because of cloud cover and/or sprinkler irrigation 14 . Farmers' adaption of the system is still limited. Additionally, the high cost of obtaining these images for relatively small areas is a significant drawback 15 . Multispectral cameras mounted on unmanned aerial vehicles (UAVs) have enormous potential to resolve this problem. UAVs can be deployed rapidly and frequently for image acquisition, resulting in reduced costs, greater flexibility in terms of data resolution and mission timing [16][17][18] . For instance, a variable rate N fertilization map was created using hyperspectral airborne images 19 , the ground-sensor measurements were compared with hyperspectral images to determine the N sufficiency index 17 , and estimation of N side-dress using NDVI computed from aerial imaging 20 . However, radiometric and geometric calibrations are needed for the UAVs on-board miniaturized electro-optical sensors to obtain quantitative results and provide precise georeferencing 21 . UAVs also fail to perform on-board mosaicking images due to limited computational resources. Laser induced breakdown spectroscopy (LIBS) is an analytical method for qualitative and quantitative elemental detection. LIBS can be readily applied to soil samples, providing rapid chemical analysis of soil samples and their constituents (e.g., Nitrogen, Potassium, Phosphorus, Calcium). The combination of an autonomous UA, LIBS, and machine learning can be used to achieve in-field measurement which provides instant results for deficient nutrient analysis and fertilization planning. With appropriate calibration, the LIBS analysis can provide quantitative measurement for most elements in soil including, carbon, nitrogen, potassium, sulfur, and phosphorus 22,23 . There have been some applications of standalone LIBS systems in precision agriculture [22][23][24][25] . However, there has been no detailed research of LIBS application in combination with ML and UAVs. Some studies have found it challenging to measure nitrogen using LIBS due to environmental factors; Earth's atmosphere is almost 80% nitrogen which will interfere with the sample measurement result since the soil is less than 1% nitrogen. Testing in a vacuum or in low-pressure conditions has been suggested to improve measurement accuracy 23 . In this study, we conducted LIBS analysis on soil samples under a normal atmosphere for observation. Low laser pulse energies were used to minimize the breakdown of air and thereby minimize the influence of atmospheric nitrogen.
The purpose of the present study is to develop a machine learning (ML)-based predictive model to estimate TN of soil using crops and soil spectral characteristics measured from the multispectral images captured from a UAV, and LIBS. Specifically, we train a multi-layer perceptron regression (MLP-R) and support vector regression (SVR) model to predict TN in soil. We use root mean square error (RMSE) and computational time (CT) as performance metrics to measure the performance of the above predictive model. To reduce the RMSE and lower CT in the machine learning models, we perform hyper-parameter optimization (HPO). The HPO tuning process depends on the ML model used for prediction 26 . The traditional way to tune hyper-parameter is through manual testing, although it requires a deep understanding of the ML models 27 . However, manual tuning is ineffective for many problems due to a large number of hyperparameters, model complexity, time-consuming model evaluations, and non-linear hyper-parameter interactions. Several HPO techniques 28 have been used for different applications such as grid search (GS), random search (RS), bayesian optimization, genetic algorithm (GA), and particle swarm optimization. In this study, we implement GS, RS, and GA for hyper-parameter optimization.

Experimental design
An aerial survey was carried out with Mavic 2 Pro UAV. We obtained multispectral images using the Sentera high-precision NDVI single sensor which was mounted on the UAV (Fig. 1a). The sensor is 1.2 MP CMOS with a 60 • horizontal FOV and a 47 • vertical FOV and works with two wide spectral bands: red (625 nm CWL × 100 nm width) and NIR (850 nm CWL × 40 nm width) with a pixel count of 1248 horizontal/950 vertical. The green band is typically unused. The sensor has a total weight of 30 g and size of 25  www.nature.com/scientificreports/ multispectral images at each of the growth stages (V4, V8, and V12) and used Sentera image stitching software to mosaicking the images. The multispectral images were captured while the UAV was following the raster scan pattern using the following parameters and experimental setup, (i) Parameters: • • Desired resolution: The Ground Sample Distance (GSD)/pixel of the multispectral camera was set to 0.05 m for a 60.96 m altitude. • Cloud cover and time of day: The UAV was flown when the sun was highest in the sky for more accurate data. Data is best when sky conditions are consistent, ideally 100% sunny or 100% cloudy. Flying with a mix of sun and clouds causes inconsistency in brightness and contrast while stitching images. Therefore, the stitched image will provide an inaccurate NDVI value.
We collected six soil samples 6.1 m from the edge of the field and six samples from the opposite side of the field, and six soil samples from the center of the field as shown in Fig. 1b. We followed the soil sampling methods for South Dakota region 29 to select the sample locations and number of samples collected. A total of 54 soil samples were collected at an 0.2 m depth from six patches (3 samples per patch) at the V4, V8, and V12 stages (18 samples per stage) using a hydraulic probe. We avoided sampling from the areas where conditions were different from the rest of the field (e.g., former manure piles, fertilizer bands, or fence lines). Figure 1b shows the sample locations across the patches.
Calibration. LIBS utilizes a high energy pulsed laser which generates a high temperature ranging from 10 • -20,000 • K resulting in plasma formation when focused on a sample. This, in turn, leads to ablation of a minuscule amount of sample, leading to excitation of the sample's constituent elements. As the plasma cools, these excited atoms and electrons emit photons which correspond to specific elements present in the sample. These photons are collected by a spectrometer and result in quantitative and qualitative analysis of samples. The SciAps Z-300 handheld LIBS analyzer was used for these measurements. This device has an extended spectrometer wavelength range from 190 µm to 950 µm . The extended range allows emission lines from elements H, F, N, O, Br, Cl, Rb and, S to be measured. The LIBS instrument is equipped with a Q-switched Nd:YAG laser, 5-6 mJ per pulse at 1064 nm. Ten laser pulses are shot on the soil samples in the presence of Ar purge to obtain averaged data on each measurement. The focused laser on the soil surface forms a µm size of a sample into > 10,000 • K plasma. The unique emission spectrum is collected by the spectrometer as the plasma cools.
We used NIST LIBS database 30 to determine the N lines from the emission spectrum (Fig. 2) and found N lines at 493.4 nm, 746.6 nm, 821.4 nm, and 868.1 nm (Fig. 3). However, we discarded the 746.6 nm N lines due to weaker intensity response and inconsistency between the samples in wavelength. We verified the N lines from the study of soil nutrient detection for precision agriculture 22 . From the soil samples, we select four samples randomly and obtained the actual TN of soil in ppm for calibration. We analyzed all the 54 soil samples in LIBS to determine the N spectrum's maximum intensity at 493.4 nm, 821.4 nm, and 868.1 nm (Figs. 5, 6, and 7) at V4, V8, and V12 stages.
Using the correlation between actual TN and the maximum intensity of N spectrum, we construct calibration plots for 493.4 nm, 821.4 nm, and 868.1 nm through linear regression (Fig. 4). We use R 2 as our calibration metric and find R 2 = 0.98 , R 2 = 0.99 , and R 2 = 0.90 , respectively, showing a strong correlation between the actual soil TN and the peak intensity of the N spectrum. Using the calibrated model, we converted the peak intensity of the N spectrum (Figs. 5, 6, and 7) to TN (ppm) for all the 54 soil samples (Table 1) to generate the training data for the ML models.   (1) and (2), band separation (Fig. 8a) was performed to compute NDVI (Fig. 8b) and extract the pixel values from each of the bands. The dataset (Table 1)    www.nature.com/scientificreports/

Methods
In supervised learning, the goal is to obtain an optimal predictive model function f * based on the input x and the output y to minimize the cost function L(f(x), y). In this study, we particularly use MLP-R and SVR which can be used for both classification and regression problems. We applied HPO techniques to determine the best set of hyper-parameters from the ML models and train the ML models using those hyper-parameters on the training dataset.   Table 2) was created using the solver type 33 , activation function 34 , learning rate and hidden layer sizes.

Support vector regression (SVR). Support vector machine (SVM) makes data points linearly separa-
ble by mapping them from low-dimensional to high-dimensional space. The classification boundary creates a partition between the data points by generating a hyperplane 35 . SVM concepts can be applied to regression problems by generalizing them. SVR uses a symmetrical loss function that penalizes both high and low misestimates equally. The ε-tube is used to generalize SVM to SVR by adding an ε-insensitive region around the function, ignoring the absolute values of errors less than a certain threshold ε from both above and below the estimation 36,37 . In SVR, points outside the tube are penalized, but points inside the tube, whether above or below the function, are not penalized. SVR uses different types of kernels for non-linear functions to map the data into a higher dimensional space 34,36,37 . Linear kernels, radial basis function (RBF), polynomial kernels, and sigmoid kernels are common kernel types in SVR. We created the hyper-parameter configuration ( Table 2) using the kernel types, regularization parameter (C) 34 , and distance error (epsilon) of the loss function 34 .
Hyper-parameter optimization (HPO). GS, RS, and GA HPO techniques were executed within their respective hyper-parameters to train the model. We performed cross-validation by splitting the train and test data into 5-folds. After obtaining the RMSE from the cross-validation score, we selected the hyper-parameters which yielded the lowest RMSE. Finally, using the best set of hyper-parameters we trained the MLP-R and SVR  www.nature.com/scientificreports/ Table 1. Generated training data for the ML models. www.nature.com/scientificreports/ models for each HPO technique. Figure 9b shows the step-by-step process of HPO, training the dataset, and prediction of test data. GS exhaustively evaluates all the combinations in the hyper-parameter configuration space specified by the user in the form of a grid configuration 38 . The user must identify the global optimums manually since GS cannot utilize the well-performing regions 28 . However, in RS, the user defines a budget (i.e., time) as well as the upper and lower bounds of the hyper-parameter values. RS randomly selects the values from the pre-defined boundary and trains until the budget is exhausted 28 . If the configuration space is wide enough, RS can detect the global optima. Assuming a model has k parameters and each of them has n distinct values, the GS computational complexity increases exponentially at a rate of O(n k ) 39 . Therefore, the effectiveness of GS depends on the size of the hyper-parameter configuration space. For RS, the computational complexity is defined as O(n), where n is specified by the user before the optimization process starts 28 .

Data count Red (DN) NIR (DN) Green (DN) NDVI RH (%) Air temp (C)
GA 40 randomly initializes the population and chromosomes. Genes represents the entire search space, hyperparameters, and hyper-parameter values. GA uses a fitness function to evaluate the performance of each individual in the current generation similarly to the objective function of a ML model. To produce a new generation, GA performs selection, crossover, and mutation operations on the chromosomes involving the next hyperparameter configurations to be evaluated. The cycle continues until the algorithm reaches the global optimum.

Results and discussion
To evaluate the HPO methods, we implemented five fold cross-validation and used RMSE as the performance metric. Additionally, we measured CT as a model efficiency metric. CT is the total time required to complete an HPO process. We specify the same hyper-parameter configuration space (Table 2) for all HPO methods to fairly compare GS, RS and GA. The optimal hyper-parameter configuration (Table 3) was determined by each of the HPO methods based on the lowest RMSE for all three wavelengths. We tuned the models on a machine with an 8 Core i7-9700K processor and 16 gigabytes (GB) of memory. We used Python 3.5, multiple open-source Python libraries, and open-source Python frameworks, including sklearn 34 . Figure 10 shows that for both MLP-R and SVR, RS produces much faster results than GS while maintaining lower RMSE for the same search space size. In general, GA offers lower RMSE for both models but has  www.nature.com/scientificreports/ a higher CT compared to GS and RS in all three wavelengths. Overall, MLP-R outperforms SVR in terms of performance. However, we achieved better efficiency with SVR in our dataset. We introduced the machine learning approach to estimate the TN of soil using NDVI and multispectral characteristics (R, NIR and G) of the images. We also consider the environmental factors such as air temperature and RH. The performance of MLP-R and SVR models were tested on a fixed configuration space for the hyperparameters under various hyper-parameter optimization techniques at three different wavelengths (Table 4). For both MLP-R and SVR, the default HP configuration do not yield the lowest RMSE, this demonstrates the significance of utilizing HPO. From Table 4, the estimation error of predicting soil TN is lowest in GA compared to GS and RS for both MLP-R and SVR, where µ is the mean and σ is the standard deviation. While training the models, we split our dataset into train and test for all three wavelengths individually, where we use 80% of the data for training and 20% for testing.
The UMS framework can be used to estimate the total nitrogen in soil. However, depending on the types of soil and crops, the model needs to be re-calibrated. More specifically, the actual TN of soil should be obtained from the subset of the samples to calibrate the N spectrum's intensity after determining the N lines using LIBS. Furthermore, N lines that fall around the 500 nm region should be avoided in sea sand due to the interferences with Titanium (Ti) lines 23 .

Conclusions
In this paper, we have demonstrated the ability of a UAV-based multispectral sensing solution to estimate soil total nitrogen. Specifically, we implemented two machine learning models multilayer perceptron regression and support vector regression to predict soil total nitrogen using a suite of data classes including UAV-based imaging data in red, near infrared, and green spectral bands, normalized difference vegetation indices (computed using the multispectral images), air temperature, and relative humidity. We performed hyperparameter optimization methods to tune the models for prediction performance. Overall, our numerical studies confirm that our machine learning-based predictive models can estimate total nitrogen of the soil with a root mean square percent error (RMSPE) of 10.8%.

Data availability
The source code, and the training data can be found here, https:// git. io/ JOaqK.