Improving Biomass and Grain Yield Prediction of Wheat Genotypes on Sodic Soil Using Integrated High-Resolution Multispectral, Hyperspectral, 3D Point Cloud, and Machine Learning Techniques

of


Introduction
Sodic soils occupy 581 million hectares globally, representing one of the significant constraints to agricultural production [1,2]. Worldwide, Australia has the widest cover of sodic soils (340 million hectares) [2]. Wheat is an important rain-fed crop grown in northeastern Australia, and sodicity poses a major threat to its production as it limits plant available water in the root zone [3][4][5]. To increase farm productivity and profitability, one important strategy is to identify wheat genotypes with better tolerance to sodic soil conditions. To do this, it is of importance to be able to evaluate and compare the performance of genotypes in terms of their growth and yield response in sodic soil [6].
Crop phenotyping, i.e., the evaluation and quantitative measurement of complex traits based on different characteristics over different growth stages [7], is widely used to identify adaptive traits and quantify plant-soil-environment interactions [8,9]. Phenotyping has traditionally relied on the collection of information about crop growth (e.g., biomass, leaf appearance, crop height, crop nutrient content) over the growing season using manual field sampling methods. However, manual methods are often time-consuming, labor-intensive, and potentially inaccurate [10,11]. The use of unmanned aerial vehicles fitted with a variety of powerful and high-resolution sensors has provided an alternative method of data collection that is inexpensive, as well as time-and labor-efficient, and which can provide accurate soil-crop information in real time at a field to regional scale [12][13][14]. However, limited studies have been reported using these approaches for crop phenotyping in constrained environments, particularly on sodic soils. These sensor-based approaches may be useful for phenotyping crops/cultivars grown on sodic soils where spatial variability and incomplete canopy cover can be major challenges to representative sampling. Although a recent study reported the potential of a high-resolution UAV-thermal imaging sensor to evaluate physiological performance, water status, and growth of wheat cultivars on sodic soils [15], greater research is required using different cropping traits and multiple sensor-based approaches for phenotyping on sodic soils that can improve estimation of crop growth and yield for adaptation of wheat genotypes in sodic soil environments.
Remote sensing spectral information and traits obtained from visible and near-infrared (VNIR) wavelengths can be used to determine crop health, vigor, moisture content, and growth [16]. With recent advances in the airborne LIDAR and/or UAV red-green-blue (RGB) sensor-based 3D point cloud techniques in high-throughput phenotyping (HTP), it is also possible to obtain high-resolution crop architectural traits (height, volume, and canopy cover) [17,18]. These architectural traits have been reported as important in influencing aboveground biomass and/or grain yield variation [19][20][21]. In addition, a wide range of proximal and handheld instruments, such as EMI (Geonics EM38 ® , Genonics LTD, Canada) [22,23], GreenSeeker ® handheld sensor [24], and ASD FieldSpec ® range [25,26], are now available to provide sound, detailed information on soil and crop characteristics at close distances. Hyperspectral sensors with hundreds of spectral bands have demonstrated potential for assessing and predicting crop performance and yield [26]. Researchers have suggested that using combined trait information derived from a number of sensors can be advantageous to quantify complex crop information on non-sodic soil [27,28]. However, the effectiveness of these remote sensing techniques has only been reported in the absence of soil constraints. These combined traits information may, thus, also be useful for assessing the performance of crops growing in sodic soils. However, there is currently little published information available to identify appropriate physiological traits that can be used to evaluate and forecast the relative performance of crops and/or cultivars using these proximal and handheld instruments with that of UAVs on sodic soils. Hence, the identification of robust techniques is required for phenotyping on sodic soils. A wide range of vegetation indices (VIs), such as normalized difference vegetation index (NDVI), enhanced vegetation index (EVI), and optimized soil adjusted vegetation index (OSAVI), can be computed from multispectral and/or hyperspectral data collected using UAVs or proximal sensing. These VIs provide information about crop health, greenness, and vigor [29][30][31]. The NDVI is one of the most well-known VIs for monitoring crop growth and predicting grain yield [32][33][34][35][36]. However, sparse vegetation or soils that generate high reflectance (especially in dryland agriculture) can adversely affect the reliability of NDVI. Rondeaux et al. [37] reported that OSAVI was an improved model of NDVI that could measure canopy reflectance by using 0.16 as an optimal value for reducing the variation of canopy reflectance in the presence of soil background conditions. In contrast, Gao et al. [38] reported that EVI can be more responsive to canopy cover, structural variation, and architecture in comparison to NDVI, especially at a regional scale [39]. Huete et al. [40] also suggested that EVI decouples the background effects from the canopy and is less sensitive to atmospheric influence. Furthermore, NDVI and EVI together can be used to detect temporal vegetation cover and provide insights into canopy biophysical information. While a number of different VIs exist to predict crop growth and yield, it is often difficult to determine which is the most useful for a specific study, and the relative performance of VIs can vary depending on the site-specific soil and environmental conditions. There is also limited information on the usefulness of VIs in crop growth and yield estimation, particularly on sodic soils. Hence, a comprehensive assessment of these VIs derived from a different remote and/or proximal sensors is required to select appropriate traits for reliable crop growth and yield forecasting on sodic soils.
An accurate and reliable seasonal crop growth and yield forecast is essential for agriculture risk assessment and improving productivity. A number of multivariate predictive machine learning (ML) algorithms can be used in agricultural yield forecasting and modeling [41][42][43][44]. ML has been proven to improve decision making and provide reliable outcomes that are consistent and robust [45][46][47]. Although a number of studies have reported the usefulness of various ML algorithms, including support vector machine (SVM), artificial neural network (ANN), classification and regression tree (CRT), Gaussian process regression (GPR), and multilinear regression (MLR), for agricultural yield estimation at various spatial extensions (Table 1), the authors are currently unaware of any published assessment on the usefulness of these integrated ML and optical remote sensing-based approaches for crop growth and yield estimations, especially in constrained sodic soils. Hence, a thorough evaluation of these ML approaches is required to determine its utility for crop growth and yield prediction on sodic soils.
Accurate crop growth and yield prediction can be reliant on the nature and importance of input cropping traits and data variability. It can also be specific to different environments and sites. Hence, a thorough evaluation of these approaches is necessary for adaptation in sodic soils. The present study proposes an integrated optical remote sensing and ML-based framework to quantify wheat genotypic traits to assist estimation of growth and yield on sodic soils. The research objectives are to (1) assess the utility of a variety of remote sensing techniques and traits to determine crop growth and yield of wheat genotypes on sodic soils, (2) examine the abilities of different ML algorithms coupled with remote sensing-derived biophysical traits to improve prediction of wheat biomass and grain yields, and (3) quantify the impacts of sodic soils on wheat crop growth and yield of to identify tolerant genotypes on sodic soils. The present study tests if integrated optical remote sensing and ML-based techniques can accurately forecast yield and help to distinguish tolerant genotypes on sodic soil that can assist agronomists/researchers to select appropriate techniques/traits for phenotyping on sodic soils by reducing the need for extensive, labor-intensive, and manual methods of phenotyping to examine the traits and crop growth. In addition, the study may also guide decisions for farmers to select genotypes tolerant to sodic soil constraints to maximize productivity in sodic soil environments.

Site Selection and Soil Sampling
The experiment was performed at two rain-fed sites, one moderately sodic (MS) (28.15 • S and 150.22 • E) and the other highly sodic (HS) (28.08 • S and 150.15 • E), near Goondiwindi in northeastern Australia. Both sites have well-structured gray vertisols with high clay content and are located at an average elevation of 268 m above mean sea level. During the normal crop growing season from May to October, air temperature varied betweeñ 5 and~35 • C for both sites with seasonal means of 14.7 • C for the MS and 14.9 • C for the HS site [15,54]. Mean relative humidity was 54.6% for the MS and 54.3% for the HS site, and total in-season crop rainfall was 86.0 mm for the MS and 85.6 mm for the HS site ( Figure 1) [15].

Site Selection and Soil Sampling
The experiment was performed at two rain-fed sites, one moderately sodic (MS) (28.15°S and 150.22°E) and the other highly sodic (HS) (28.08°S and 150.15°E), near Goondiwindi in northeastern Australia. Both sites have well-structured gray vertisols with high clay content and are located at an average elevation of 268 m above mean sea level. During the normal crop growing season from May to October, air temperature varied between ~5 and ~35 °C for both sites with seasonal means of 14.7 °C for the MS and 14.9 °C for the HS site [15,54]. Mean relative humidity was 54.6% for the MS and 54.3% for the HS site, and total in-season crop rainfall was 86.0 mm for the MS and 85.6 mm for the HS site ( Figure 1) [15]. In-season crop rainfall and air temperature (monthly mean) of the moderately sodic (MS) and highly sodic (HS) sites in the wheat-growing season of 2018 [15]. The error bars represent the standard error of the monthly mean values.
At each site, soil samples were collected from a minimum of eight points using a hydraulic sampling rig to take 50 mm diameter soil cores to a depth of 150 cm. The soil samples were dried at 40 °C and ground to <2 mm. In a 1:5 soil water suspension [55], pH, EC, and chloride (Cl) were measured; the EC of saturated extract (ECse) was then computed from EC 1:5, and the clay content was determined using the pipette method [56,57]. Exchangeable Na and the cation exchange capacity (CEC) were measured using a 1 M NH4Cl (pH 8.5) extraction solution [58], and ESP was calculated from exchangeable Na relative to CEC. The volumetric moisture content as a percentage was determined after drying samples at 105 °C to obtain the gravimetric soil moisture content before multiplying this by the soil bulk density (BD) [15]. The soil physicochemical properties of both sites are presented in Figure 2. In-season crop rainfall and air temperature (monthly mean) of the moderately sodic (MS) and highly sodic (HS) sites in the wheat-growing season of 2018 [15]. The error bars represent the standard error of the monthly mean values.
At each site, soil samples were collected from a minimum of eight points using a hydraulic sampling rig to take 50 mm diameter soil cores to a depth of 150 cm. The soil samples were dried at 40 • C and ground to <2 mm. In a 1:5 soil water suspension [55], pH, EC, and chloride (Cl) were measured; the EC of saturated extract (EC se ) was then computed from EC 1:5, and the clay content was determined using the pipette method [56,57]. Exchangeable Na and the cation exchange capacity (CEC) were measured using a 1 M NH 4 Cl (pH 8.5) extraction solution [58], and ESP was calculated from exchangeable Na relative to CEC. The volumetric moisture content as a percentage was determined after drying samples at 105 • C to obtain the gravimetric soil moisture content before multiplying this by the soil bulk density (BD) [15]. The soil physicochemical properties of both sites are presented in Figure 2.

Experimental Design and Crop Biophysical Measurements
We used a randomized complete block design (RCBD) for the experiment with eight replications [15,54,59]. Four of the replications (72 plots) were used for destructive plant sampling, crop height, and biophysical measurements close to flowering stage, 110-112 days after sowing (DAS), and another four replications (72 plots) were used for grain yield measurements at maturity (152 DAS). Eighteen contrasting wheat genotypes were tested (Supplementary Table S1) [15,54,59], making a total of 144 plots (four columns and 36 rows) ( Figure 3). Destructive sampling plots and grain harvest plots for each genotype in each replicate block were arranged in adjacent pairs. Each plot was 5 × 2 m at harvest and consisted of five planting rows with 30 cm spacing in between (Supplementary Figure S1) [15,54,59]. The crops were sown in late May and harvested in early November of 2018. The aboveground biomass yield at both sites was harvested from a 3 × 0.5 m area from the middle three rows of each destructive sampling plots at 110-112 DAS (Supplementary Figure S1), stored in paper bags, and then dried at 70 °C for 72 h to determine plant dry weight (DW) [15,48]. In addition, crop heights were measured manually by placing a ruler vertically on the soil surface and measuring the height of plants from the middle three rows in each plot. Average plant height of a plot was then estimated using the mean of these measurements. Grain yield was measured from the yield plots by machine

Experimental Design and Crop Biophysical Measurements
We used a randomized complete block design (RCBD) for the experiment with eight replications [15,54,59]. Four of the replications (72 plots) were used for destructive plant sampling, crop height, and biophysical measurements close to flowering stage, 110-112 days after sowing (DAS), and another four replications (72 plots) were used for grain yield measurements at maturity (152 DAS). Eighteen contrasting wheat genotypes were tested (Supplementary Table S1) [15,54,59], making a total of 144 plots (four columns and 36 rows) ( Figure 3). Destructive sampling plots and grain harvest plots for each genotype in each replicate block were arranged in adjacent pairs. Each plot was 5 × 2 m at harvest and consisted of five planting rows with 30 cm spacing in between (Supplementary Figure S1) [15,54,59]. The crops were sown in late May and harvested in early November of 2018.

Experimental Design and Crop Biophysical Measurements
We used a randomized complete block design (RCBD) for the experiment with eight replications [15,54,59]. Four of the replications (72 plots) were used for destructive plant sampling, crop height, and biophysical measurements close to flowering stage, 110-112 days after sowing (DAS), and another four replications (72 plots) were used for grain yield measurements at maturity (152 DAS). Eighteen contrasting wheat genotypes were tested (Supplementary Table S1) [15,54,59], making a total of 144 plots (four columns and 36 rows) ( Figure 3). Destructive sampling plots and grain harvest plots for each genotype in each replicate block were arranged in adjacent pairs. Each plot was 5 × 2 m at harvest and consisted of five planting rows with 30 cm spacing in between (Supplementary Figure S1) [15,54,59]. The crops were sown in late May and harvested in early November of 2018. The aboveground biomass yield at both sites was harvested from a 3 × 0.5 m area from the middle three rows of each destructive sampling plots at 110-112 DAS (Supplementary Figure S1), stored in paper bags, and then dried at 70 °C for 72 h to determine plant dry weight (DW) [15,48]. In addition, crop heights were measured manually by placing a ruler vertically on the soil surface and measuring the height of plants from the middle three rows in each plot. Average plant height of a plot was then estimated using the mean of these measurements. Grain yield was measured from the yield plots by machine harvesting at crop maturity (152 DAS). The aboveground biomass yield at both sites was harvested from a 3 × 0.5 m area from the middle three rows of each destructive sampling plots at 110-112 DAS (Supplementary Figure S1), stored in paper bags, and then dried at 70 • C for 72 h to determine plant dry weight (DW) [15,48]. In addition, crop heights were measured manually by placing a ruler vertically on the soil surface and measuring the height of plants from the middle three rows in each plot. Average plant height of a plot was then estimated using the mean of these We used a handheld GreenSeeker ® 505 (Trimble Inc., Sunnyvale, CA, USA) active sensor to measure NDVI from tillering to crop maturity to acquire data on growth and greenness during these crop development stages. NDVI values range between 0 and +1. Higher values indicate greater canopy cover and/or greenness suggesting good canopy health. The sensor was conveyed at an average speed of 0.84 m/s and an elevation of 0.5 m from the top of the canopy from the middle three rows of each biomass sampling plot (Supplementary Figure S1). The measurements were taken continuously and at the same speed without stopping the sensor in between plots [15]. The field of view of the sensor was oval, and it changed with the elevation from the ground. Bare soil reflectance was measured separately from outside of the plot area and calibrated using a calibration panel to reduce the effect of soil reflectance on the canopy [15,59].
Canopy reflectance was also measured using a spectroradiometer (ASD FieldSpec ® HandHeld 2, Malvern Panalytical Ltd., Malvern, UK) under cloud-free, wind-free, and sunny conditions between 11:00 a.m. and 3:00 p.m. Weekly in-situ monitoring of crop development was conducted, and a date near flowering (110-112 DAS) was selected for crop assessment. A date near flowering was chosen since this is the most crucial time before grain filling when the canopy is most fully developed and likely to modify the reflectance most efficiently [15]. The instrument was calibrated using a standard white spectralon calibration panel at each time of point shoot to reduce the effects of soil reflectance on canopy spectral measurements. The continuous spectra were recorded in Vis to NIR from 325 nm to 1075 nm with ±1 nm accuracy and a resolution of <3 nm for each spectrum. The field of view of the sensor was 25 • and, in each plot, five spectral measurements of 0.2 m diameter area of the canopy were recorded from the middle three rows (Supplementary Figure S1) by locating the sensor in a vertical position 0.5 m above the canopy cover [26]. The measured canopy reflectance data were then exported into an ASCII file using ViewSpec Pro, an integrated package with RS 3 spectral acquisition software (Analytical Spectral Devices Inc., Boulder, CO, USA) [60]. Subsequently, the raw data were exported into the R statistical software platform and converted into hyperspectral data by using the hyperSpec package [61]. The preprocessing ('cleaning') was carried out using a spreadsheet program. Accordingly, the bands introducing excessive variation ('noise') into the data (901-1050 nm) and very short wavelength bands (325-399 nm) were excluded.

UAV-Based Sensing
The multispectral data were obtained at both sites via UAV campaigns ( Figure 4) at 110-112 DAS in sunny, cloud-free, and low-wind conditions during midday. The propeller aero points were positioned on the ground uniformly for recording GPS positioning information ( Figure 4d) before starting the flight. The flight plan was set up in the Pix4D capture software platform (Pix4D S.A., Switzerland) using a 2D polygon mission tool with lower ground sample distance (GSD = 1.3 cm/px), 85% image overlapping, a field of view (FOV) of 85 • , and 30 m flight height. A five-band multispectral RedEdge-M camera (MicaSense Inc., Seattle, DC, USA) ( Figure 4c) with a 5.4 mm focal length was used to capture high-resolution (1280 × 960 pixels) multispectral imagery. The camera was mounted on a DJI Matrices 100 (SZ DJI Technology, Shenzhen, Guangdong, China) quadcopter. The imagery was acquired with bandwidth and the center wavelength of 20 nm at 475 nm (blue), 20 nm at 560 nm (green), 10 nm at 668 nm (red), 10 nm at 717 nm (red edge), and 40 nm at 840 nm (near infrared) [34,62]. The raw images were recorded in 16 bit raw GeoTiff files stored on a local digital card, with the images uploaded to ATLAS cloud (MicaSense Inc., Seattle, DC, USA) for visualization, further processing, and information extraction. soil reflectance using Otsu's thresholding method [64,65] in the MATLAB R2020a image processing platform (The Mathworks® Inc., Natick, MA, USA) (Supplementary Figure  S2), which finds the optimal threshold by maximizing the weighted sum of between-class variances. Lastly, statistics of plot-wise mean canopy reflectance were extracted from the canopy pixels using the zonal statistics tool in Arc GIS 10.7.1 software platform (ESRI Ltd., CA, USA). In addition, a high-resolution RGB camera (DJI Technology, Shenzhen, Guangdong, China), fitted with a 24 mm lens with a field of view (FOV) of 25° × 20° and 5472 × 3648 pixel resolution, was used to acquire RGB images of the site at 110-112 DAS. The camera was mounted on a DJI Phantom 4 Professional UAV at 30 m flight height with a GSD of 1.3 cm/px, 85% image overlapping, and a FOV of 85°. The RGB mission planning and setup was similar to that used for the UAV multispectral mission. The images were then loaded and aligned in Agisoft Metashape 1.5.5 software platform (Agisoft LLC, St. Petersburg, Russia) for orthomosaicing and further preprocessing and point cloud extraction. Initially, a sparse point cloud was generated with the formation of camera positioning [66]. Gradually, dense point clouds were generated using a dense stereo-matched algorithm with the structure from motion (SfM) technique and exported for crop height extraction in the Lidar 360 software platform (GreenValley International, USA). The average point density of the SfM captures was 824 and 945 pt/m 2 for the MS and the HS site, respectively. Furthermore, an * .LAS file (Lidar data exchange) was created and converted to a multipoint structure for measuring crop height. Lastly, by taking a horizontal cross-sectional profile of the point cloud for each plot ( Figure 5), the vertical difference ('z' direction) between the reference soil surface and the apex point (average top of the terminal spikelet) was measured [67], and plot-wise crop heights were calculated by taking the mean of at least the 5000 highest cloud data points in each plot (Supplementary Figure S3). Image preprocessing, including orientation correction, tie point matching, reflectance calibration, and orthomosaicing was performed in Agisoft Metashape 1.5.5 software platform (Agisoft LLC, St. Petersburg, Russia) [63]. A Calibrated Reflectance Panel (CRP) (MicaSense Inc., Seattle, DC, USA) was used before and after the individual flights for calibrating GeoTiff layers of five spectral bands for the conversion of digital number (DN) into reflectance (Figure 4e). Pixels of the canopy were segmented from background soil reflectance using Otsu's thresholding method [64,65] in the MATLAB R2020a image processing platform (The Mathworks ® Inc., Natick, MA, USA) (Supplementary Figure S2), which finds the optimal threshold by maximizing the weighted sum of between-class variances. Lastly, statistics of plot-wise mean canopy reflectance were extracted from the canopy pixels using the zonal statistics tool in Arc GIS 10.7.1 software platform (ESRI Ltd., Redlands, CA, USA).
In addition, a high-resolution RGB camera (DJI Technology, Shenzhen, Guangdong, China), fitted with a 24 mm lens with a field of view (FOV) of 25 • × 20 • and 5472 × 3648 pixel resolution, was used to acquire RGB images of the site at 110-112 DAS. The camera was mounted on a DJI Phantom 4 Professional UAV at 30 m flight height with a GSD of 1.3 cm/px, 85% image overlapping, and a FOV of 85 • . The RGB mission planning and setup was similar to that used for the UAV multispectral mission. The images were then loaded and aligned in Agisoft Metashape 1.5.5 software platform (Agisoft LLC, St. Petersburg, Russia) for orthomosaicing and further preprocessing and point cloud extraction. Initially, a sparse point cloud was generated with the formation of camera positioning [66]. Gradually, dense point clouds were generated using a dense stereomatched algorithm with the structure from motion (SfM) technique and exported for crop height extraction in the Lidar 360 software platform (GreenValley International, Berkeley, CA, USA). The average point density of the SfM captures was 824 and 945 pt/m 2 for the MS and the HS site, respectively. Furthermore, an *.LAS file (Lidar data exchange) was created and converted to a multipoint structure for measuring crop height. Lastly, by taking a horizontal cross-sectional profile of the point cloud for each plot (Figure 5), the vertical difference ('z' direction) between the reference soil surface and the apex point (average top of the terminal spikelet) was measured [67], and plot-wise crop heights were calculated by taking the mean of at least the 5000 highest cloud data points in each plot (Supplementary Figure S3).

Vegetation Indices Derived from UAV Multispectral and Proximal Hyperspectral Data
Several vegetation indices (NDVI, OSAVI, and EVI) were derived from UAV multispectral imagery and proximal hyperspectral data to evaluate crop characteristics ( Table  2). The NDVI was calculated as the simple normalization process between near-infrared (NIR) and red (R) spectral bands using Equation (1). We also computed OSAVI as an improved VI over the NDVI. The OSAVI enhances canopy reflectance by subduing the soil reflectance using an optimized soil adjusted coefficient (L= 0.16) in Equation (2). The EVI was computed using Equation (3) for its ability to optimize crop reflectance with improved sensitivity for the biomass yield prediction.

Regression Analysis
A linear regression model was fitted in MATLAB R2020a (The Mathworks® Inc., Natick, MA, USA) to predict biomass and grain yield in response to individual VIs and crop heights using a 10-fold cross-validation method. A test of significance at p < 0.05 for a difference from '0' was used to determine significance of biomass and grain yield prediction. The accuracy of the models was assessed by measuring the root-mean-square error (RMSE) (Equation (4)) and coefficient of determination (R 2 ) Equation (5) in the cross-validation. In general, small RMSE and high R 2 values can be used to determine the accuracy of a model.

Vegetation Indices Derived from UAV Multispectral and Proximal Hyperspectral Data
Several vegetation indices (NDVI, OSAVI, and EVI) were derived from UAV multispectral imagery and proximal hyperspectral data to evaluate crop characteristics ( Table 2). The NDVI was calculated as the simple normalization process between near-infrared (NIR) and red (R) spectral bands using Equation (1). We also computed OSAVI as an improved VI over the NDVI. The OSAVI enhances canopy reflectance by subduing the soil reflectance using an optimized soil adjusted coefficient (L = 0.16) in Equation (2). The EVI was computed using Equation (3) for its ability to optimize crop reflectance with improved sensitivity for the biomass yield prediction. Table 2. Vegetation indices evaluated in this study. [40,70,71] NIR is near infrared, R is red, B is blue wavelengths, G is the gain factor = 2.5, C 1 = 6 and C 2 =7.5 are the coefficients of aerosol resistance term, and L = 1 is the canopy background adjustment factor in Equation (3). For hyperspectral VIs, we used the reflectance of the wavebands 864 nm (NIR), 671 nm (R), and 467 nm (B) [40,68,70].

Regression Analysis
A linear regression model was fitted in MATLAB R2020a (The Mathworks ® Inc., Natick, MA, USA) to predict biomass and grain yield in response to individual VIs and crop heights using a 10-fold cross-validation method. A test of significance at p < 0.05 for a difference from '0' was used to determine significance of biomass and grain yield prediction. The accuracy of the models was assessed by measuring the root-mean-square error (RMSE) (Equation (4)) and coefficient of determination (R 2 ) Equation (5) in the crossvalidation. In general, small RMSE and high R 2 values can be used to determine the accuracy of a model.
where n is the total number of data points, z i represents the observed value at i, z(x) i is the predicted value at i, z i represents an average observed value, and z(x) i represents an average predicted value.

Analysis of Variance (ANOVA)
We performed two-way ANOVA using the SPSS ® Statistics 25 software platform (IBM, New York, NY, USA) to determine whether any significant differences existed among the site means of wheat genotypes based on VIs, 3D point cloud-derived crop height, biomass yield (g/m 2 ), and grain yield (g/m 2 ) at the MS and HS sites. The testing of significant differences between levels within each factor was performed using Fisher's protected least significant difference (LSD) test at the 5% significance level. The variables were treated as nested effects, and replicates were fitted as a random effect.

Machine Learning
We used four multivariate supervised regression learning techniques in ML, i.e., MLR, SVM, GPR, and ANN, to predict biomass yield and grain yield for both MS and HS sites as a function of site-specific VIs and 3D point cloud-derived crop heights in MATLAB R2020a environment (The Mathworks ® Inc., Natick, MA, USA). These ML algorithms were selected specifically, since previous studies on non-sodic soils (see Table 1) have identified the diversified usefulness of these techniques in agricultural applications, including yield prediction, and they can also be applied on small to moderate datasets, for example, at a small field scale. A principal component analysis (PCA) was used for model feature selection to understand explained data variance in the models before applying the ML algorithms for prediction. The ML architecture employed in this study is illustrated in Figure 6 with a stepwise process including the design of training and testing datasets, feature selection, cross-validation, model performance evaluation, and decision making.

Multitarget Linear Regression
We used and compared linear, robust linear, and stepwise algorithms under this category. In general, a multi-robust linear regression algorithm is less sensitive to outliers and also provides flexibility in fitting [72]. The stepwise regression uses variables (VIs, crop height, biomass yield, and grain yield) at each step and removes variables that are not significant (p > 0.05) for the prediction [73].

Support Vector Machine Regression
We used and compared different kernel functions in SVM, such as linear, Gaussian, cubic, and quadratic for their ability to accurately predict the outcomes. The SVM is considered a reliable supervised learning technique that is used for classification, pattern recognition, regression, and prediction [74]. The linear kernel is used for a linearly separable dataset [75], whereas the Gaussian SVM is used for complex relationships [76]. Quadratic and cubic functions can provide flexibility in fitting [72].

Gaussian Process Regression
The GPR can predict the outcomes accurately. However, it is often difficult to interpret the outcomes. As recommended by The MathWorks Inc. [72], the preset flexibility in model type gallery was maintained in the MATLAB R2020a environment to avoid overfitting and to produce a small error in training. The kernel functions, such as exponential, squared exponential, matern 5/2, and rational quadratic were used and compared in their ability to predict the outcomes. Remote Sens. 2020, 17, x FOR PEER REVIEW 11 of 28

Gaussian Process Regression
The GPR can predict the outcomes accurately. However, it is often difficult to interpret the outcomes. As recommended by The MathWorks Inc. [72], the preset flexibility in model type gallery was maintained in the MATLAB R2020a environment to avoid overfitting and to produce a small error in training. The kernel functions, such as exponential, squared exponential, matern 5/2, and rational quadratic were used and compared in their ability to predict the outcomes.

Artificial Neural Network
A multilayer perceptron model (MLP) is the most used network model in ANN. The backpropagation algorithm of MLP is used to search for the lowest error function in the weight space using gradient-based methods [77]. To avoid overfitting in the MLP result, it is important to train the network model with an appropriate number of hidden layers and neurons [78]. In this network model, a sum of inputs and bias passes through a transfer function is used to deliver the output to the activation point. The perceptron network is composed of artificial nodes or neurons. These nodes or neurons are the units, which can process information in layers and are connected by synaptic weights. They can transfer information in a supervised way for building a prediction model after the classification of data. The classified data are stored in a memory. Typically, the model is interconnected by the nodes or neurons in a three-layer network structure (input, hidden, and output layer). Each node or neuron of a layer is linked to the next layer. The input data carry information on the network structure, which is multiplied by weights while entering the hidden layer. Furthermore, the hidden layer processes this information, and the output layer classifies and predicts an outcome. The weights are the preassigned numbers, which are in the form of an argument that is passed to a nonlinear mathematical function and an activation function to return numbers between 0 and 1 [79].
A multilayer feedforward network model was performed in this study for prediction of biomass and for prediction of grain yield (Figure 7). Each multilayer feedforward network model consisted of a number of input, hidden, and output layers. The network model was selected by optimizing hidden layers, neurons, and functions using multiple runs and cross-validated by assigning 70% of the dataset to training, 15% to testing, and the remaining 15% to model validation. The network was trained using a Levenberg-Marquardt backpropagation algorithm. This algorithm is known to be used for its quick convergence in backpropagation [80,81]. The perceptron determines the prediction outputs by evaluating the input weights.
fer function is used to deliver the output to the activation point. The perceptron n is composed of artificial nodes or neurons. These nodes or neurons are the units can process information in layers and are connected by synaptic weights. They ca fer information in a supervised way for building a prediction model after the classi of data. The classified data are stored in a memory. Typically, the model is interco by the nodes or neurons in a three-layer network structure (input, hidden, and layer). Each node or neuron of a layer is linked to the next layer. The input da information on the network structure, which is multiplied by weights while ente hidden layer. Furthermore, the hidden layer processes this information, and the layer classifies and predicts an outcome. The weights are the preassigned numbers are in the form of an argument that is passed to a nonlinear mathematical function activation function to return numbers between 0 and 1 [79].
A multilayer feedforward network model was performed in this study for pr of biomass and for prediction of grain yield (Figure 7). Each multilayer feedforw work model consisted of a number of input, hidden, and output layers. The n model was selected by optimizing hidden layers, neurons, and functions using m runs and cross-validated by assigning 70% of the dataset to training, 15% to testi the remaining 15% to model validation. The network was trained using a Levenbe quardt backpropagation algorithm. This algorithm is known to be used for its qu vergence in backpropagation [80,81]. The perceptron determines the prediction by evaluating the input weights. The performance of the MLs was evaluated by measuring RMSE (Equation (Equation (5)), mean squared error (MSE) (Equation (6)), and mean absolute error (Equation (7)) in the cross-validation. As mentioned earlier in Section 2.5.1, a smal MSE, and MAE for biomass and grain yield prediction was used to determine t mum performance of the ML models.
where n is the total number of data points, represents the observed value at , the predicted value at , ̅ represents an average observed value, and ̅ repre average predicted value. The performance of the MLs was evaluated by measuring RMSE (Equation (4)), R 2 (Equation (5)), mean squared error (MSE) (Equation (6)), and mean absolute error (MAE) (Equation (7)) in the cross-validation. As mentioned earlier in Section 2.5.1, a small RMSE, MSE, and MAE for biomass and grain yield prediction was used to determine the optimum performance of the ML models.
where n is the total number of data points, z i represents the observed value at i, z(x) i is the predicted value at i, z i represents an average observed value, and z(x) i represents an average predicted value.

Soil Constraints and Agro-Climatic Conditions
The depth-wise distribution of soil physicochemical properties ( Figure 2) showed that both sites have very similar clay contents (41-58%) (Figure 2e) and, thus, are likely to have a similar water holding capacity. Both sites also have a higher ESP in the subsoil (18.2% for the MS and 24% for the HS site, at 100-150 cm) than the surface soil (2.7% for the MS and 14% for the HS site, at 0-10 cm) (Figure 2d). In addition, both sites have high Cl concentrations in the subsoil at 100-150 cm, with the HS site (>2300 mg/kg) having substantially higher Cl than the MS site (700-750 mg/kg) (Figure 2b).
Initial (prior to sowing) volumetric soil moisture data (Figure 2f) showed that both the sites had adequate and mostly similar soil moisture content stored at 0-150 cm depth (a mean of~36.6% for the MS and~38.8% for the HS site), which helped in the germination of crops despite low in-crop rainfall during sowing to emergence (~10 mm for both MS and HS sites) (Figure 1). Agro-climatic conditions at both sites were similar (p > 0.01) from May to October 2018 (Figure 1), suggesting that crops at both sites received similar rainfall during the growing season.

GreenSeeker ® NDVI
Crop stagewise in situ monitoring from the measurement of GreenSeeker ® NDVI during crop development indicated that maximum canopy development occurred at approximately 110-112 DAS (just before flowering) at both sites (Supplementary Figure S4). At this time, the NDVI for healthy cultivars reached a peak of~0.74 for the MS and~0.53 for the HS site; thus, this was identified as a suitable time to differentiate between cultivars.
Higher NDVI values at the MS site also indicated more favorable crop growth and canopy greenness compared to the HS site (Supplementary Figure S4). Previous studies reported that, in non-constrained soils, the NDVI values for healthy vegetation gradually increase from the beginning of the crop season and reach a peak at the middle of the season or before flowering [82,83]. The data from this study also suggest that plants growing on a constrained soil follow a similar pattern of development, with close to flowering (110-112 DAS) being the most suitable time for obtaining remote sensing data of crops to evaluate and differentiate crop conditions for forecasting biomass and grain yields. Furthermore, biomass and grain yield were both significantly and positively correlated with GreenSeeker ® NDVI near flowering at both sites (Table 3). Table 3. Cross-validation of parameters for biomass yield and grain yield as a function of GreenSeeker ® -measured NDVI near flowering on a moderately sodic (MS) and a highly sodic (HS) site. Coefficient of determination (R 2 ) values are presented with different asterisks to depict probabilities of a significant difference from '0' with a significant level of p < 0.05 *, p < 0.001 **, and p < 0.0001 ***; n = 72.

UAV Multispectral Imaging
The VIs derived from UAV multispectral data were also significantly and positively correlated with biomass and grain yield at both sites ( Figure 8). When comparing the UAV multispectral VIs, the EVI was more closely associated with biomass yield at both sites (R 2 = 0.82; RMSE = 41.6 g/m 2 for the MS and R 2 = 0.67; RMSE = 40.0 g/m 2 for the HS site) compared to the OSAVI and NDVI (Figure 8 and Table 4), whereas the NDVI was more closely correlated with grain yield at both sites (R 2 = 0.75; RMSE = 17.7 g/m 2 for the MS and R 2 = 0.53, RMSE = 25.4 g/m 2 for the HS site) compared to the EVI and OSAVI (Figure 8 and Table 4).

Proximal Hyperspectral Sensing
Non-imaging and ground-based proximal hyperspectral sensing-derived VIs showed a significant positive association with biomass and grain yield at both sites ( Figure 9). Similar to the UAV multispectral VIs, the hyperspectral data-derived EVI was more closely associated with biomass yield at both sites (R 2 = 0.56; RMSE = 65.4 g/m 2 for the MS and R 2 = 0.42; RMSE = 53.3 g/m 2 for the HS site) compared to the OSAVI and NDVI (Figure 9 and Table 4), whereas the NDVI had a closer agreement with grain yield at both sites (R 2 = 0.48; RMSE = 25.5 g/m 2 for the MS and R 2 = 0.37; RMSE = 29.2 g/m 2 for the HS site) than the EVI and OSAVI (Figure 9 and Table 4).

Figure 8.
Relationship between UAV-multispectral imagery-derived VIs with biomass yield and grain yield at a moderately sodic (MS) (a,b) and highly sodic (HS) (c,d) site; the matrix shows the coefficient of determination (R 2 ) values between the variables, as also indicated by the color scale (right). Coefficient of determination (R 2 ) values are presented with different asterisks to depict probabilities of a significant difference from '0' with a significant level of p < 0.001 ** and p < 0.0001 ***; n = 72.

Proximal Hyperspectral Sensing
Non-imaging and ground-based proximal hyperspectral sensing-derived VIs showed a significant positive association with biomass and grain yield at both sites (Figure 9). Similar to the UAV multispectral VIs, the hyperspectral data-derived EVI was more closely associated with biomass yield at both sites (R 2 = 0.56; RMSE = 65.4 g/m 2 for the MS and R 2 = 0.42; RMSE = 53.3 g/m 2 for the HS site) compared to the OSAVI and NDVI ( Figure  9 and Table 4), whereas the NDVI had a closer agreement with grain yield at both sites (R 2 = 0.48; RMSE = 25.5 g/m 2 for the MS and R 2 = 0.37; RMSE = 29.2 g/m 2 for the HS site) than the EVI and OSAVI (Figure 9 and Table 4). The VIs derived from handheld, proximal hyperspectral sensing exhibited a range of generally higher values compared to UAV multispectral VIs, indicating a slightly greater variability in handheld hyperspectral data within a site compared to those using UAV. For example, the maximum peak of the proximal hyperspectral NDVI value reached~0.92 and~0.71 at the MS and HS site, while the UAV-derived multispectral NDVI reached the maximum peak of~0.72 and~0.49 at the MS and HS site (110-112 DAS). This might be due to the changes in prevailing environmental conditions in the field during a longer period of survey using a handheld, proximal sensor that might have caused larger data variability compared to the rapid and automated UAV-based sensing [84,85]. Furthermore, UAV multispectral VIs were more closely associated with both biomass and grain yield than ground-based proximal GreenSeeker ® NDVI or ASD FieldSpec hyperspectral VIs, as indicated by the R 2 values.
Remote Sens. 2020, 17, x FOR PEER REVIEW 15 of 28 Figure 9. Relationship between non-imaging, ground-based proximal hyperspectral data-derived VIs with biomass yield and grain yield at a moderately sodic (MS) (a,b) and highly sodic (HS) (c,d) site; the matrix shows the coefficient of determination (R 2 ) values between the variables, as also indicated by the color scale (right). Coefficient of determination (R 2 ) values are presented with different asterisks to depict probabilities of a significant difference from '0' with a significant level of p < 0.05 *, p < 0.001 **, and p < 0.0001 ***; n = 72.  . Relationship between non-imaging, ground-based proximal hyperspectral data-derived VIs with biomass yield and grain yield at a moderately sodic (MS) (a,b) and highly sodic (HS) (c,d) site; the matrix shows the coefficient of determination (R 2 ) values between the variables, as also indicated by the color scale (right). Coefficient of determination (R 2 ) values are presented with different asterisks to depict probabilities of a significant difference from '0' with a significant level of p < 0.05 *, p < 0.001 **, and p < 0.0001 ***; n = 72.

UAV RGB Sensor-Based 3D Point Cloud Techniques
Biomass and grain yield were significantly and positively correlated with 3D point cloud-derived crop height, which explained 71% and 50% of the variability of biomass yield for the MS and HS site and 56% and 39% of the variability in grain yield for the MS and HS sites, respectively (Table 5), which was a slightly better association compared to manual, ground-based measurements of crop height at both sites (Table 5). On the other hand, a significant association (p < 0.01) was observed between ground-measured crop height and 3D point cloud crop height for both sites (R 2 = 0.53 and 0.44 for the MS and HS site, respectively) (Supplementary Figure S5); a lower association between ground-measured crop height and biomass and/or grain yield suggests a larger data variability in groundmeasured crop height compared to 3D point cloud crop height. This is sensible, since ground-measured crop height using a ruler relies on manual eye estimation of the field observer, whereas high-throughput UAV measurements, automated data processing, and extraction of crop height using software algorithms might have improved data quality and precision. This suggests that there may be a reduced need for extensive, labor-intensive, and manual measurements of crop height in the field. There is little or no suggestion that, for indicating biomass and yield potential, manual height measurements are superior to data obtained by UAV. Table 5. Cross-validation of model parameters for biomass yield and grain yield as a function of crop heights at a moderately sodic (MS) and a highly sodic (HS) site. Coefficient of determination (R 2 ) values with different asterisks indicate different probabilities that the correlation differed from '0' with a significant level of p < 0.05 *, p < 0.001 **, and p < 0.0001 ***; n = 72.

Comparing ML Algorithms for Prediction of Biomass and Grain Yield on Rain-Fed Sodic Soil
We compared different ML algorithms and different kernel functions integrated with UAV multispectral VIs and 3D point cloud crop height to test if these improve the estimation of biomass and grain yield of wheat on rain-fed sodic soil. The smallest prediction error likely suggests the best-performing model.
Under multitarget regression, the stepwise kernel was found to be effective in the prediction of both biomass and grain yield and have the lowest error estimates compared to the multiple linear and/or multi-robust linear kernels (Tables 6 and 7). In SVM, the linear kernel function achieved superior prediction accuracy (lowest error) in comparison to quadratic, cubic, and Gaussian kernel functions at both MS and HS sites (Tables 6 and 7). In GPR, the squared sequential, matern 5/2, and rational quadratic functions predicted biomass yield with slightly less error than an exponential function. However, the exponential kernel achieved superior accuracy in grain yield prediction in comparison to the other GPR kernels (Tables 5 and 6). Overall, the MLP kernel in ANN achieved greatest estimated accuracy in the prediction of both biomass and grain yields at the two sites compared to all ML algorithms and kernels used in this study. The VIs and crop height combinedly explained 89% and 82% of the variability of biomass yield at the MS and HS site and 88% and 74% of the variability of grain yield at the MS and HS site using an MLP kernel in ANN (Tables 6 and 7). The best-fitted models using ANN are illustrated in Supplementary  Figures S6 and S7. The results suggest that canopy spectral information and crop height are useful indicators of aboveground biomass and/or grain yield of wheat on rain-fed sodic soil, and they can improve our understanding and ability to forecast biomass yield and grain yield variability in-field. Table 6. Performance evaluation between ML algorithms and kernels based on biomass yield prediction as a function of combined UAV multispectral VIs and 3D point cloud crop height on sodic soil; n = 72.

Model feature selection
PCA explained a total of 95% of variance.

Comparing Crop Growth, Biomass, and Grain Yields on Sodic Soils
The ANOVA outcome showed that the VIs and 3D point cloud-derived crop height was significantly different (p < 0.001) at MS and HS sites. The mean NDVI, OSAVI, EVI, and crop height were significantly higher at the MS site than at the HS site ( Figure 10). As they are important indicators of growth and health of crops, significantly lower values at the HS site suggest that crops grown at the HS site experienced more stress than at the MS site. This stress was primarily due to sodic soil constraints rather than environmental factors, as the in-season crop rainfall, air temperature, and initial moisture stored in the soil profile were similar at both sites (Figures 1 and 2).
When comparing the performance of wheat genotype biomass yield between sites, the current study indicated that most of the genotypes had significantly higher predicted biomass (p < 0.05) at the MS site than the HS site close to flowering (Figure 11), whereas, at maturity, genotypic grain yield performance was substantially lower (p < 0.01) at the HS site compared to the MS site ( Figure 12). The results suggest that crop stress most likely increased with development and had a significant negative influence on grain yield at the maturity stage, especially for the genotypes at the HS site.
When considering a single site, the discrimination between genotypes was more prominent at the HS site than the MS site (Figure 12), suggesting variable effects of high sodic soil constraints on tolerance and performance of wheat genotypes. This also indicates that the techniques used in this study can be useful to differentiate between cultivars in more sodic soil conditions. The study identified that wheat genotypes Mitch, Corack, Mace, Trojan, Lancer, and Bremer were more tolerant to sodic soil constraints than Emu Rock, Janz, Flanker, and Gladius. However, seasonal changes and changes in soil constraint level may have diversified influences over genotypic performance. Hence, implementation of these strategies over different seasonal and sodic soil conditions, using various of crops and/or cultivars, is recommended for future studies. Overall, the results support the suggestion that integrated UAV optical remote sensing and ML techniques have the ability to improve estimation of crop growth and yield on sodic soils that can be used to identify cultivars tolerant to sodic soil environments.
Remote Sens. 2020, 17, x FOR PEER REVIEW 19 of 28 level may have diversified influences over genotypic performance. Hence, implementation of these strategies over different seasonal and sodic soil conditions, using various of crops and/or cultivars, is recommended for future studies. Overall, the results support the suggestion that integrated UAV optical remote sensing and ML techniques have the ability to improve estimation of crop growth and yield on sodic soils that can be used to identify cultivars tolerant to sodic soil environments.   level may have diversified influences over genotypic performance. Hence, implementation of these strategies over different seasonal and sodic soil conditions, using various of crops and/or cultivars, is recommended for future studies. Overall, the results support the suggestion that integrated UAV optical remote sensing and ML techniques have the ability to improve estimation of crop growth and yield on sodic soils that can be used to identify cultivars tolerant to sodic soil environments.

Discussion
The present study aimed to develop an integrated optical remote sensing and MLbased framework to quantify wheat genotypic traits to assist the estimation of growth and yield on sodic soils. We tested if integrated optical remote sensing and ML-based techniques can accurately forecast yield and help to distinguish tolerant genotypes on sodic soil, which can assist agronomists/researchers to select appropriate techniques/traits for phenotyping on sodic soils by reducing the need for extensive, labor-intensive, and manual methods of phenotyping to examine the traits and crop growth.

Traits and Sensor Performances
Fern et al. [86] reported that the OSAVI can better estimate green biomass and vegetative cover than NDVI in non-constrained semi-arid rangeland soil. Similarly, our study indicated that both EVI and OSAVI are slightly better indicators of biomass yield than NDVI in rain-fed sodic soil environment (Figures 8 and 9 and Table 4). However, a number of studies reported that NDVI can be a useful indicator of grain yield performance on non-sodic soils [35,[87][88][89]. Likewise, the present research found that NDVI is a more useful indicator of wheat grain yield than OSAVI and EVI in the presence of sodic soil constraints (Figures 8 and 9 and Table 4). Crop rainfall and ambient air temperature can have a large influence on canopy vigor and greenness, particularly in rain-fed conditions. Adequate rainfall between late tillering and flowering crop growth stages (Figure 1) in the current study could have allowed us to obtain sufficient canopy spectral cues that were closely correlated with biomass and grain yield.
When comparing sensors, the high-resolution UAV RedEdge-M multispectral sensor showed a slightly greater performance in estimating crop growth and yield than groundbased proximal hyperspectral and GreenSeeker ® sensors (Figures 8 and 9). A systematic and rapid UAV flight plan and operation might have potentially helped to achieve a good quality of imaging data with the multispectral sensor when compared to handheld proximal sensors, as UAV multispectral imagery-derived traits were more closely associated with both biomass and grain yields than either proximal GreenSeeker or ASD FieldSpec sensor-derived VIs (Table 3 and Figures 8 and 9). Saura et al. [90] also reported that a high

Discussion
The present study aimed to develop an integrated optical remote sensing and MLbased framework to quantify wheat genotypic traits to assist the estimation of growth and yield on sodic soils. We tested if integrated optical remote sensing and ML-based techniques can accurately forecast yield and help to distinguish tolerant genotypes on sodic soil, which can assist agronomists/researchers to select appropriate techniques/traits for phenotyping on sodic soils by reducing the need for extensive, labor-intensive, and manual methods of phenotyping to examine the traits and crop growth.

Traits and Sensor Performances
Fern et al. [86] reported that the OSAVI can better estimate green biomass and vegetative cover than NDVI in non-constrained semi-arid rangeland soil. Similarly, our study indicated that both EVI and OSAVI are slightly better indicators of biomass yield than NDVI in rain-fed sodic soil environment (Figures 8 and 9 and Table 4). However, a number of studies reported that NDVI can be a useful indicator of grain yield performance on non-sodic soils [35,[87][88][89]. Likewise, the present research found that NDVI is a more useful indicator of wheat grain yield than OSAVI and EVI in the presence of sodic soil constraints (Figures 8 and 9 and Table 4). Crop rainfall and ambient air temperature can have a large influence on canopy vigor and greenness, particularly in rain-fed conditions. Adequate rainfall between late tillering and flowering crop growth stages (Figure 1) in the current study could have allowed us to obtain sufficient canopy spectral cues that were closely correlated with biomass and grain yield.
When comparing sensors, the high-resolution UAV RedEdge-M multispectral sensor showed a slightly greater performance in estimating crop growth and yield than groundbased proximal hyperspectral and GreenSeeker ® sensors (Figures 8 and 9). A systematic and rapid UAV flight plan and operation might have potentially helped to achieve a good quality of imaging data with the multispectral sensor when compared to handheld proximal sensors, as UAV multispectral imagery-derived traits were more closely associated with both biomass and grain yields than either proximal GreenSeeker or ASD FieldSpec sensorderived VIs (Table 3 and Figures 8 and 9). Saura et al. [90] also reported that a high precision of UAV data can be achieved with a systematic flight plan in automatic mode, and they demonstrated the potential benefits of using it for improved crop yield prediction.
The accuracy of the handheld spectroradiometer used to acquire the hyperspectral data might have been reduced due to weather and prevailing environmental conditions. For example, the longer time taken to obtain handheld measurements compared to the flight time of the drone would have allowed the environment to change more during the period of measurement. The spectral data combine a mixture of diffuse and specular reflectance because of different particle sizes, scattering, and multicollinearity [91]. Variation in spectral reflectance is often caused by differences in the size of particles in the object. This causes a deviation of light at various angles as a function of the wavelengths and creates scattering effects [84]. Ana Belén, Enoc, Víctor Marcelo, Marta, and José Ramón [84] also reported that results did not satisfactorily improve even after normalization and preprocessing of hyperspectral spectroradiometer data for a vineyard study. However, our study was unable to pinpoint the exact reason for variations in the data. The problem might be overcome by using UAV hyperspectral sensing with an automated flight plan, as UAV sensing can capture information rapidly and, thus, reduce the influence of environmental variability. Pádua et al. [85] observed that unmanned aerial systems (UAS) with specific imaging sensors achieve better results than handheld devices.
Yang and Jinfei [92] estimated canopy height using UAV 3D point cloud data and demonstrated the ease with which the technique can be used to obtain crop growth information at the field level. A few studies reported that UAV-based crop height measurement may be integrated with cropping traits to improve the estimation of the aboveground biomass of barley and/or winter wheat on non-constrained soil [53,93,94]. Hence further research is warranted to employ these integrated sensor-based techniques and traits to improve estimation of biomass and/or yield on rain-fed sodic soils. UAV RGB 3D point cloud techniques show promise for the extraction of crop height, which is closely related to aboveground biomass and grain yield on sodic soils.

Yield Prediction on Rain-Fed Sodic Soils Using ML
We compared four ML algorithms (MLR, SVM, GPR, and ANN) for crop growth and yield prediction on sodic soils and observed that the ANN predicted aboveground biomass and grain yield with slightly less error than the others. The results suggest that the ANN has superior ability to account for site-specific complexities of agricultural traits that can be used to improve estimation of growth and yield, particularly on sodic soils. Recently, ANN has become popular for different applications in agriculture, including crop yield prediction. For example, Safa et al. [95] reported a better ability of an ANN model relative to MLR for estimating wheat yield using a heterogenous dataset in multiple arable farms on non-sodic soils in New Zealand. Mokarram and Bijanzadeh [96] also reported a better prediction accuracy using the MLP (ANN) network model in comparison to MLR for barley biomass yield (R 2 = 0.89) and grain yield (R 2 = 0.92) on non-constrained soils. Recently, Dhakal, Gautam, and Bhattarai [76] reported that an ANN model with 10 neurons in the hidden layer performed comparatively better in training with R 2 = 0.84 and RMSE = 1.5 MJ·m −2 ·day −1 when estimating global solar radiation. Furthermore, their study also found that the stepwise kernel has good potential for estimation with less error (R 2 = 0.88 and RMSE = 1.5 MJ·m −2 ·day −1 ). Ekici et al. [97] reported that the Matern 5/2 kernel in GPR produced better results in comparison to other supervised ML algorithms and kernels for a hybrid power system study. However, there are a limited number of studies on agricultural biomass and grain yield prediction using UAV-and ground-based optical sensor-derived VIs, using crop height as input crop parameters with comprehensive ML approaches, particularly on sodic soil. The performance of the models can be influenced by sample size, site characteristics, areal coverage, and the traits used for the growth and yield prediction study. On sodic soils, the main challenge lies in accurate plant sampling and extraction of useful cropping traits due to larger spatial variability and incomplete crop cover caused by soil constraints compared to non-constrained soils [48].
Hence, the appropriate selection of traits and the accurate extraction of crop information using high-throughput sensing and techniques were essential to feed the ML models for improved prediction accuracy. It was also necessary to test and compare model performance individually at the MS and HS site with varying levels of soil constraints to investigate the performance of complex traits for evaluating crop growth and yield variation. Crop production on sodic soil has previously been reported to be negatively affected due to high Cl concentration in subsoil, which decreased water and nutrient uptake by root [5,98]. However, those were labor-intensive and time-consuming techniques to determine the relative performance of crops/cultivars and tolerance to soil constraints [5,98]. Hence, the present study demonstrated the efficacy of optical remote sensing and statistical modeling approaches to help identify appropriate traits. Overall, our results support the suggestion that integrated optical UAV remote sensing and ML techniques have the potential to improve estimation of crop growth and yield on sodic soils. This information can also be important in the discrimination between the performance of cultivars/genotypes.

Crop Growth and Yield Vary with Changes in Sodic Soil Constraints
Recent studies have suggested that constraints in sodic soil can increase crop water stress and negatively affect the physiological performance of wheat [15,48]. A study reported that the UAV thermal imaging-derived crop water stress index had a close negative agreement with wheat growth and yield at critical crop development stages, such as close to flowering and maturity on rain-fed sodic soils [48], suggesting that wheat grown on sodic soils appears water-stressed. In this study, we also observed a significant variation in spectral indices and crop growth close to flowering across different sodic soil conditions ( Figure 10) and the variation was most likely due to water limitations. High subsoil ESP and Cl concentrations have a strong potential to restrict root water and nutrient uptake from the soil and reduce rooting depth [98,99]. Thus, in both soils, crops would be expected to be stressed and have reduced growth, with profound adverse impacts at the HS site, which supports the suggestions made by an earlier study that high sodic soil constraints can reduce wheat physiological growth [15]. Studies have also reported that canopy temperature is an important trait to evaluate wheat physiological growth on sodic soils under water-limited environments [15,48]. Results from the current study indicated that optical remote sensing-based spectral VIs and crop height information are also strong indicators of wheat biomass and grain yields on rain-fed sodic soils and can be used to differentiate crop growth, stress, and yield performance between different levels of sodicity. Overall, the results confirm the adverse impacts of soil constraints, which can have a large influence on biomass and/or grain yield (Figures 11 and 12). Genotypic biomass and yield performance variations/differences were observed within a site, and it was clearer between the sites with the changes in soil constraint level. The results closely correspond to the findings of a recently published study on sodic soil [48], where the authors demonstrated that water stress can further increase with the seasonal development of wheat crop in highly sodic soil environments compared to moderately sodic soils and can adversely affect wheat yield at maturity. Overall, the results support the suggestion that wheat genotypic growth and yield performance vary with sodic soil constraints, and that genotypic tolerance to soil constraints can be determined using the appropriate selection of traits and integrated high-throughput approaches, including UAV optical remote sensing and ML.

Conclusions
This study successfully evaluated the potential of optical remote sensing and ML algorithms to improve the estimation of growth and yield of wheat genotypes, which is useful to differentiate genotypic tolerance to sodic soils. A range of approaches were tested and proposed to demonstrate their application to crops grown in sodic soil. The potential of a number of sensors and ML algorithms was assessed for their ability to accurately forecast wheat biomass and grain yields for identifying constraint tolerant cultivars on sodic soils. Different VIs derived from a UAV multispectral (RedEdge-M) sensor predicted biomass yield and grain yield at least as well as, and possibly better than, the handheld proximal hyperspectral spectroradiometer and/or GreenSeeker ® sensors. Although hyperspectral sensing was reported to be superior in previous studies conducted on non-constrained soils, we could not achieve high accuracy using a handheld non-imaging hyperspectral sensor. Crop height extracted from UAV RGB-based 3D point cloud was found to have a closer association with both biomass and grain yield than ground-based manually measured crop height, suggesting the adaptation of this technique on sodic soil. Furthermore, ML algorithms integrated with optical remote sensing-derived traits were tested and compared for their ability to improve forecasting of crop growth and yield on sodic conditions. Overall, the present study demonstrated the efficacy of optical remote sensing and ML-based techniques in providing real-time and detailed information on appropriate trait selection for evaluating crop growth and yield variation to identify constraint-tolerant cultivars in a rain-fed sodic soil environment. Follow-up studies should focus on implementing these techniques with a greater diversity of genotypes and at multiple sodic sites to compare the results. This study can greatly assist agronomists/plant physiologists in the selection of appropriate techniques/traits for phenotyping on sodic soils by reducing the need for extensive, labor-intensive, and manual methods of phenotyping to examine the relative performance of crop traits, growth, and yield. In addition, this may help guide decisions by farmers and breeders in selecting tolerant varieties in sodic soil environments to sustain agricultural production.
The key research findings are the following: • High sodic soil constraints negatively affected crop growth and development and reduced yield. • A number of the methods were able to discriminate differences between sites and some between genotypes within a site.

•
The UAV multispectral (RedEdge-M) sensor performed with slightly less error than the ground-based handheld proximal hyperspectral and/or GreenSeeker ® sensors for the measurements of crop traits.

•
The UAV RGB-based 3D point cloud technique is promising for crop height measurements and suggests there is a reduced need for the manual, labor-intensive, and tedious process of crop height measurements in the field.

•
The EVI was in more close association with biomass yield and the NDVI with grain yield on sodic soils. • Integrated VIs and crop height were useful indicators of biomass and grain yield performance of wheat genotypes on rain-fed sodic soil.

•
The ANN performed slightly better than multitarget regression, SVM, and GPR in estimating biomass yield and grain yield on sodic soils.

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/10 .3390/rs13173482/s1, Figure S1: Drone image of a representative plot illustrating a typical destructive biomass sampling area ('yellow' square box indicating 0.5 × 3 middle rows area) and plant reflectance measurements area using GreenSeeker and ASD FieldSpec ('orange' rectangular box including 3 middle rows) within the plot. One of the propeller markers for GPS location calibration is shown in the lower left hand of the image, Figure S2: A representative site image showing classified and separated soil and canopy pixels using Otsu's algorithm, Figure S3: Institutional Review Board Statement: "Not applicable" for studies not involving human or animals.