Sugarcane Nitrogen Concentration and Irrigation Level Prediction Based on UAV Multispectral Imagery

Sugarcane is the main industrial crop for sugar production, and its growth status is closely related to fertilizer, water, and light input. Unmanned aerial vehicle (UAV)-based multispectral imagery is widely used for high-throughput phenotyping, since it can rapidly predict crop vigor at field scale. This study focused on the potential of drone multispectral images in predicting canopy nitrogen concentration (CNC) and irrigation levels for sugarcane. An experiment was carried out in a sugarcane field with three irrigation levels and five fertilizer levels. Multispectral images at an altitude of 40 m were acquired during the elongating stage. Partial least square (PLS), backpropagation neural network (BPNN), and extreme learning machine (ELM) were adopted to establish CNC prediction models based on various combinations of band reflectance and vegetation indices. The simple ratio pigment index (SRPI), normalized pigment chlorophyll index (NPCI), and normalized green-blue difference index (NGBDI) were selected as model inputs due to their higher grey relational degree with the CNC and lower correlation between one another. The PLS model based on the five-band reflectance and the three vegetation indices achieved the best accuracy (Rv = 0.79, RMSEv = 0.11). Support vector machine (SVM) and BPNN were then used to classify the irrigation levels based on five spectral features which had high correlations with irrigation levels. SVM reached a higher accuracy of 80.6%. The results of this study demonstrated that high resolution multispectral images could provide effective information for CNC prediction and water irrigation level recognition for sugarcane crop.


Introduction
As the most important sugar crop, sugarcane is mainly grown in tropical and subtropical areas and provides approximately 80% of the world's sugar [1,2]. The growth of sugarcane is closely related to fertilizer, water, and radiation intensity. Evaluating the growth situation of sugarcane in a timely manner, and adjusting the field management strategy accordingly, is of great significance to the yield and quality of sugarcane. In recent years, remote sensing with spectral images at different scales has been considered an effective high-throughput phenotyping solution for predicting the growth and yield of crops.
Large-scale spectral imagery can cover an area ranging from 25 to 3600 km 2 per image. Spatial resolutions generally range from more than 1 m to tens of meters, and some data can reach 0.3 m [3,4]. They are mainly obtained by satellites, such as Sentinel [5], Gaofen historical yield records [37]. Khaki et al. proposed a convolutional neural network model called YieldNet to predict corn and soybean yield based on MODIS products [38]. Yang et al. tried to use one-year hyperspectral imagery to train a CNN classification model to estimate corn grain yield [39]. Prodhan et al. monitored drought over South Asia using a deep learning approach with 16 years of remote sensing data [40]. It can be seen that a large volume of image data, as well as ground truth data in years, were commonly needed to provide a sufficient dataset to train a deep learning network. To collect this large number of data samples is very challenging. Therefore, the dataset of deep learning is difficult to produce in some circumstances. By contrast, the traditional machine learning methods which are generally based on statistics are suitable for most of the modeling problems when relatively small number of samples are available [41,42]. Partial least squares (PLS), extreme learning machines (ELMs), backpropagation neural networks (BPNNs), support vector machine (SVM), and others, have been widely used in crop nutrient predictions. For example, Li et al. [43] used PLS to establish 12 models of fruits and seeds for rapid analysis and quality assessment. Kira et al. established a model for estimating the chlorophyll and carotenoid contents of three tree varieties based on BPNN [44]. Chen et al. constructed a BPNN model to invert rice pigment content with several spectral parameters as input [45]. Pal et al. used an ELM algorithm to classify land covers with multispectral and hyperspectral data [46]; it achieved a better classification accuracy than models established with BPNN and support vector machine (SVM), with far less computational complexity. Different machine learning methods can suit for different cases depending on variable quantity, sample quantity, and the potential relationship between inputs and output.
In this study, in order to monitor the growth status of sugarcane canopies by a high throughput method, high-resolution multispectral images of an experimental sugarcane field were obtained by a low-altitude UAV. The objectives of this study were (1) to determine the sensitive spectral features for the predictions of the CNC and irrigation levels; (2) to establish the prediction models of the CNC based on different machine learning algorithms such as PLS, BPNN, and ELM; (3) to establish classification methods of irrigation levels based on SVM and BPNN.

Study Area
The sugarcane experimental field was in Nanning, Guangxi Autonomous Region, China (latitude 22.84 • N, longitude 108.33 • E), as shown in Figure 1. From the captured multispectral image (displayed in RGB) in the right of Figure 1, it can be seen that the experimental field had 12 plots with concrete partitions. Three irrigation treatments and five fertilization treatments were applied in the field. Urea, calcium magnesium phosphate, and potassium chloride were chosen as N, phosphorus (P), and potassium (K) fertilizers, respectively. Eight plots with different irrigations and fertilizers and four blank plots without fertilizer and irrigation (denoted by BL) were set in the field. Concrete partitions at a depth of 1.2 m were built between each plot to prevent water and fertilizer infiltration. The planted seedlings were limited to 975,000 plants per hectare. The two irrigation treatments included 180 m 3 /ha (denoted by W0.6) and 300 m 3 /ha (denoted by W1.0), while the four fertilizer treatments included F1.0 (250 kg/ha of N, 150 kg/ha of P 2 O 5 , 200 kg/ha of K 2 O), F0.9 (90% of the amount of F1.0), F1.1 (110% of F1.0) and F1.2 (120% of F1.0). Water and fertilizer were applied via drip irrigation pipes. Micronutrient fertilizers were equally applied to all the plots except the blank plots. The eight plots had the same size of 20 m × 6 m, with their different treatments denoted by W0.6F0.9, W0.6F1.0, W0.6F1.1, W0.6F1.2, W1.0F0.9, W1.0F1.0, W1.0F1.1 and W1.0F1.2.
The seed canes were planted on 24 March 2018. The seedling fertilizers, which accounted for 30% of the total fertilizer application, were applied on 11 May (28 days after planting). The tillering fertilizers, which accounted for 70% of the total fertilizer application, were applied on 29 June (67 days after planting). The irrigation schedule is listed in Table 1. Rainfall was another way that water entered the open field. The rainfall in this field was 509.8 mm from the day of planting (24 March) to the day of image acquisition (11 July), and the monthly average rainfall was 127 mm. The meteorological conditions of the experiment field, including precipitation (without the irrigation), temperature and mean relative humidity, are shown in Figure 2. It can be seen that there was almost no rainfall for 15 days before the day of image acquisition. As such, the last event of water input in a large amount was the controlled irrigation on 4 July, which was a week before canopy image acquisition. This means that the rainfall had a very limited influence on the remote evaluation of water stress conditions under specific irrigation amounts. The study site and the field management layout with different irrigation and fertilization levels based on the captured multispectral image (displayed in RGB). W0.6 represents the irrigation rate of 180 m 3 /ha, W1.0 represents the irrigation of 300 m 3 /ha; F1.0 represents the standard fertilization rate, F0.9 represents 90% of the amount of F1.0, F1.1 represents 110% of F1.0 and F1.2 represents 120% of F1.0; BL1, BL2, BL3, BL4 indicate that four blank plots without fertilizer and irrigation.
The seed canes were planted on 24 March 2018. The seedling fertilizers, which accounted for 30% of the total fertilizer application, were applied on 11 May (28 days after planting). The tillering fertilizers, which accounted for 70% of the total fertilizer application, were applied on 29 June (67 days after planting). The irrigation schedule is listed in Table 1. Rainfall was another way that water entered the open field. The rainfall in this field was 509.8 mm from the day of planting (24 March) to the day of image acquisition (11 July), and the monthly average rainfall was 127 mm. The meteorological conditions of the experiment field, including precipitation (without the irrigation), temperature and mean relative humidity, are shown in Figure 2. It can be seen that there was almost no rainfall for 15 days before the day of image acquisition. As such, the last event of water input in a large amount was the controlled irrigation on 4 July, which was a week before canopy image acquisition. This means that the rainfall had a very limited influence on the remote evaluation of water stress conditions under specific irrigation amounts.

Data Collection
The multispectral images were captured at noon on 11 July 2018 (109 days after planting), in the elongating stage. The weather was sunny, cloudless, and windless. The image acquisition system was mainly composed of a drone modeled Phantom 4 Pro (DJI, Shenzhen, China) and a multispectral image sensor RedEdge-MX (MicaSense, Seattle, WA, USA), as shown in Figure 3a,b, respectively. RedEdge-MX image sensor has five spectral bands at 475 nm (blue, B), 560 nm (green, G), 668 nm (red, R), 717 nm (red edge, RE), and 840 nm (near infrared, NIR), and is equipped with a light intensity sensor and a reflectance correction panel (Group VIII, USA, Figure 3c) for radiation correction. The optical intensity sensor can correct the influence caused by changes in sunlight on the spectral images during a flight, and the fixed reflectance correction panel can be used for reflectance transformation. The drone flew at an altitude of 40 m, with 85% forward overlap and 85% side overlap. The time interval of image acquisition was 2 s, and the ground sample distance (GSD) was 2.667 cm. Four calibration tarps with reflectivity of 5%, 20%, 40% and 60%, respectively, were also placed at the open space next to the field before image acquisition, as shown in Figure 3d. Two hundred and sixty multispectral images were finally collected.

Data Collection
The multispectral images were captured at noon on 11 July 2018 (109 days after planting), in the elongating stage. The weather was sunny, cloudless, and windless. The image acquisition system was mainly composed of a drone modeled Phantom 4 Pro (DJI, Shenzhen, China) and a multispectral image sensor RedEdge-MX (MicaSense, Seattle, WA, USA), as shown in Figure 3a,b, respectively. RedEdge-MX image sensor has five spectral bands at 475 nm (blue, B), 560 nm (green, G), 668 nm (red, R), 717 nm (red edge, RE), and 840 nm (near infrared, NIR), and is equipped with a light intensity sensor and a reflectance correction panel (Group VIII, USA, Figure 3c) for radiation correction. The optical intensity sensor can correct the influence caused by changes in sunlight on the spectral images during a flight, and the fixed reflectance correction panel can be used for reflectance transformation. The drone flew at an altitude of 40 m, with 85% forward overlap and 85% side overlap. The time interval of image acquisition was 2 s, and the ground sample distance (GSD) was 2.667 cm. Four calibration tarps with reflectivity of 5%, 20%, 40% and 60%, respectively, were also placed at the open space next to the field before image acquisition, as shown in Figure 3d. Two hundred and sixty multispectral images were finally collected.

Ground Sampling and CNC Determination
Each plot was divided into three sampling areas. Each sampling area was divided into nine grids, and one plant was randomly selected to collect the first fully unfolded leaf for each grid. Nine leaves were collected to form a leaf sample for each sampling area. A total of 36 samples were finally collected, and these were immediately brought back to the laboratory for N determination. All the samples were oven-dried at 105 °C for 30

Ground Sampling and CNC Determination
Each plot was divided into three sampling areas. Each sampling area was divided into nine grids, and one plant was randomly selected to collect the first fully unfolded leaf for each grid. Nine leaves were collected to form a leaf sample for each sampling area. A total of 36 samples were finally collected, and these were immediately brought back to the laboratory for N determination. All the samples were oven-dried at 105 • C for 30 min and afterward at 75 • C for about 24 h until at a constant weight. The dried leaves were ground and weighed to 0.3 g, and the Kjeldahl method [47] was used to determine the total nitrogen (TN, %) content. The TN of those first-leaf samples, which were then considered as the CNC (%), could be calculated by Equation (1).
where % represents the unit of TN and CNC; V 1 is the consumption volume of the acid standard solution, mL; V 0 is the titration blank volume, mL; C is the concentration of the acid standard solution, mol/L; 0.014 is the 1 mol standard titration solution equivalent to the weight of N, g; m is the weight of the sample, g.

Multispectral Image Preprocessing
Pix4DMapper software (Pix4D, Prilly, Switzerland) was used to generate the mosaic image from the 260 original multispectral images, as shown in Figure 4. The mosaic image was then imported and processed in ENVI software (L3Harris Technologies, Melbourne, FL, USA). Two preprocessing steps were conducted in ENVI, including radiation correction and geometric correction.
where % represents the unit of TN and CNC; V1 is the consumption volume of the acid standard solution, mL; V0 is the titration blank volume, mL; C is the concentration of the acid standard solution, mol/L; 0.014 is the 1 mol standard titration solution equivalent to the weight of N, g; m is the weight of the sample, g.

Multispectral Image Preprocessing
Pix4DMapper software (Pix4D, Prilly, Switzerland) was used to generate the mosaic image from the 260 original multispectral images, as shown in Figure 4. The mosaic image was then imported and processed in ENVI software (L3Harris Technologies, Melbourne, FL, USA). Two preprocessing steps were conducted in ENVI, including radiation correction and geometric correction.
Radiation calibration was implemented using the radiometric correction module in ENVI. The "empirical line" method was selected, since four calibration tarps with known reflectivity were captured in the image. An empirical line was fitted by comparing the DN values and the reflectivity of the tarps. Subsequently, all the DN values in the mosaic image were able to be converted into reflectivity.
Geometric correction was conducted to eliminate the distortion. Four ground control points were selected at four corners of the field, as marked in the false-color image in Figure 4. The "image to map" function was selected to implement geometric correction with the coordinate information of the ground control points. "Nearest neighbor", which avoids introducing new pixel values, was used to resample the image to the same coordinate system (UTM projection, WGS-84 datum) as that of the ground control points. In order to extract the region of interest (ROI) out of the background, a classification method, decision tree (DT), was used to extract the sugarcane canopy from the soil, weeds, shadow, concrete, and other interfering background features. Figure 5 shows the NDVI image of the extracted canopy, and the white dots in the figure represented 36 sampling areas. To enhance sample quantity, each area was further divided into nine grids, which were approximately 1.5 m × 2.0 m in size. The average value of each grid was calculated as the spectral sample. Therefore, a total of 324 spectral samples were extracted.  Radiation calibration was implemented using the radiometric correction module in ENVI. The "empirical line" method was selected, since four calibration tarps with known reflectivity were captured in the image. An empirical line was fitted by comparing the DN values and the reflectivity of the tarps. Subsequently, all the DN values in the mosaic image were able to be converted into reflectivity.
Geometric correction was conducted to eliminate the distortion. Four ground control points were selected at four corners of the field, as marked in the false-color image in Figure 4. The "image to map" function was selected to implement geometric correction with the coordinate information of the ground control points. "Nearest neighbor", which avoids introducing new pixel values, was used to resample the image to the same coordinate system (UTM projection, WGS-84 datum) as that of the ground control points.
In order to extract the region of interest (ROI) out of the background, a classification method, decision tree (DT), was used to extract the sugarcane canopy from the soil, weeds, shadow, concrete, and other interfering background features. Figure 5 shows the NDVI image of the extracted canopy, and the white dots in the figure represented 36 sampling areas. To enhance sample quantity, each area was further divided into nine grids, which were approximately 1.5 m × 2.0 m in size. The average value of each grid was calculated as the spectral sample. Therefore, a total of 324 spectral samples were extracted.

Extraction of VIs
VIs have been widely used to qualitatively and quantitatively evaluate vegetation cover varieties and crop vigor. NDVI is the most commonly used VI, and it is also one of the important parameters closely related to crop chlorophyll and N concentration. Besides NDVI, nine other commonly used VIs (as shown in Table 2) were also selected to compare their effects on predicting the CNC. The optimal VI or a combination of VIs was used to build the prediction models of the CNC and the irrigation levels.

Grey Relational Analysis
Grey relational analysis (GRA), also called grey incidence analysis (GIA), is an important part of grey system theory, which was developed by Julong Deng [57]. At its core, it works to determine the primary and secondary relationships between various factors by calculating grey relational degree (GRD). The higher the GRD value of any two factors, the more consistent the change between those two factors. Therefore, it can be used to select the factor with the greatest influence [58]. Let the reference sequence be = ( ), = 1,2, ⋯ , and the comparison sequence be = ( ), = 1,2, ⋯ , . The GRD value between X0 and Xi is calculated by Equations (2)

Extraction of VIs
VIs have been widely used to qualitatively and quantitatively evaluate vegetation cover varieties and crop vigor. NDVI is the most commonly used VI, and it is also one of the important parameters closely related to crop chlorophyll and N concentration. Besides NDVI, nine other commonly used VIs (as shown in Table 2) were also selected to compare their effects on predicting the CNC. The optimal VI or a combination of VIs was used to build the prediction models of the CNC and the irrigation levels.

Grey Relational Analysis
Grey relational analysis (GRA), also called grey incidence analysis (GIA), is an important part of grey system theory, which was developed by Julong Deng [57]. At its core, it works to determine the primary and secondary relationships between various factors by calculating grey relational degree (GRD). The higher the GRD value of any two factors, the more consistent the change between those two factors. Therefore, it can be used to select the factor with the greatest influence [58]. Let the reference sequence be X 0 = {x 0 (k), k = 1, 2, · · · , n} and the comparison sequence be X i = {x i (k), k = 1, 2, · · · , n}. The GRD value between X 0 and X i is calculated by Equations (2) and (3).
where ρ is the identification coefficient, and its value range is 0-1, taken here as 0.5. The GRA was conducted between all the spectral features and the CNC, which were all normalized. The GRD, which was higher than 0.8, reflected that the VI had a very strong influence on the CNC.

Correlation Analysis
Correlation coefficient (R) [59] can reflect the degree of the linear correlation between two datasets. It can be calculated by Equation (4).
where n is sample size, x i and y i are the individual sample points indexed with i; x and y are the means of x i and y i for n samples.
The higher the absolute value of R, the higher the linear correlation between the two factors. It is generally considered that 0.7 ≤ |R| < 1 indicates a very high correlation, when 0.4 ≤ |R| < 0.7, it indicates a significant correlation, and when |R| < 0.4, it indicates a low correlation. Correlation analysis can be applied for multiple purposes in modeling, including: (1) to analyze the correlations between input variables and predictors to determine sensitive variables; (2) to analyze the correlations between multiple variables, during which only the variables with significant correlations should be utilized in order to simplify the complexity of the model; (3) to analyze the correlation between the predicted values of a model and the measured values, and to evaluate the effect of the model.
In this study, the correlations between the spectral features were analyzed to pick proper variables with less redundant information for CNC modeling and irrigation level classification.

Modeling Algorithms
At present, there are many machine learning algorithms. Based on previous researches, four algorithms were selected after comprehensive consideration, as shown in Table 3. PLS, BPNN and ELM were selected for CNC modeling, and the simple validation method, hold-out [60], was selected for model validation. All 324 samples were divided into calibration set and validation set according to the ratio of 7:3.
SVM and BPNN were selected for irrigation level classification. Three-fold cross validation [61] was used to produce more validation samples, and, therefore, to generate a comprehensive confusion matrix of the classification results.
The PLS algorithm builds a model by minimizing the sum of the squares of the errors. It combines the advantages of multiple linear regression, canonical correlation analysis and principal component analysis. BPNN has the characteristics of self-learning and self-adaptation, showing a strong ability to fit nonlinear functions; it also has a strong antiinterference ability and may be suitable for complex field environments. The ELM algorithm allows for the random generation of the weights and thresholds between the input layer and hidden layers; users only need to denote the number of hidden layer neurons in the whole training process. Compared with the traditional classification algorithms, ELM has a fast-learning speed and strong generalization capability. These three algorithms have different characteristics and might achieve better prediction results under different conditions or scenarios, so all three algorithms were adopted and compared for the CNC prediction in this study. The number of principal components of the PLS was 6. The training epoch, the learning rate, and the number of hidden layers of the BPNN model were 1000, 0.05, and 22, respectively. The transfer function and the number of hidden layers of the ELM model were sigmoidal function and 50, respectively.
SVM is a classic machine learning method for classification. It maps data from a lowdimensional space to a high-dimensional space through a kernel function and separates the classes with a decision surface that maximizes the margin between the classes. Thus, SVM was selected for the irrigation levels classification in this study. Due to its strong ability to perform nonlinear mapping, BPNN is suitable for not only solving fitting problems, but also classification problems, so BPNN was also selected here for comparison with SVM in the classification of irrigation levels. The penalty factor and the kernel function of the SVM model were 10 and 0.167, respectively. The training epoch, the learning rate, and the number of hidden layers of the BPNN model were 1000, 0.1, and 10, respectively.

Accuracy Assessment Metrics
R and root mean square error (RMSE) were used to evaluate the accuracies of the CNC prediction models. R was introduced in Section 2.5, and here the correlation between the predicted values and the actual values were calculated to evaluate the accuracies of the prediction models. RMSE, which was calculated by Equation (5), can directly reflect the errors of the prediction models.
where, y i andŷ i represent the estimated value and actual value for sample i, respectively. The confusion matrix is also known as the probability matrix or error matrix [62]. It is a specific matrix for visualizing algorithm performance, and is often used to evaluate the classification results. The rows in the matrix represent the actual irrigation levels and the columns represent the predicted irrigation levels. The confusion matrix is named because it can easily indicate whether multiple classes are confused (that is, one class is expected to be another class). Common indicators including producer's accuracy (PA), user's accuracy (UA), and overall accuracy (OA), can be calculated in the confusion matrix. PA refers to the ratio of the correctly classified sample numbers in a class to the actual total numbers of that class, also called true positive rate (TPR). UA refers to the ratio of the correctly classified sample numbers in a class to the classified total numbers of that class, also called positive predictive value (PPV). OA refers to the ratio of all correctly classified sample numbers to all sample numbers of all the classes. The calculation formulas of the two indicators are shown in Equations (6)- (8).

Relation Analysis between the Spectral Features and the CNC
To eliminate the influence caused by different magnitudes of each variable, the datasets were firstly normalized before GRA. In the GRA results listed in Table 4, all the VI had a GRD larger than 0.5 with the CNC, while the SRPI had the strongest grey relation of 0.94. Correlation analysis was also conducted for these ten VIs, and the results were shown in Table 5. Most of the VIs had a very high correlation (R > 0.9) with each other, except NGBDI and RVI 2 . Two variables with a higher correlation means more redundant information is contained in them, which means that selecting both as inputs could be avoided. Based on the results in Tables 4 and 5, we can find that though most of the VIs (SRPI, NPCI, RVI, MSRI, NDVI, SIPI, OSAVI and SAVI) had very high GRD (>0.8) with the CNC, they also had very high correlations (R > 0.9) with each other. Two basic rules, as follows, should be considered in the selection of spectral features, which could help by selecting the more sensitive and less-redundant spectral features as the input for CNC prediction based on the GRA and R results: (1) The GRD of the selected spectral feature(s) should be relatively high with the CNC (GRD > 0.65); (2) The R between the spectral feature(s) should be relatively low to avoid introducing redundant information and variable coupling. Therefore, not all the VIs with higher GRD could be selected as model inputs. Here, the first two VIs in Table 4, SRPI and NPCI, were recommended as efficient inputs. Furthermore, NGBDI, which had relatively higher grey relational degree with the CNC and a lower R with SRPI and NPCI, was also recommended as another efficient input variable.

Modeling with the Five-Band Reflectance
The five-band reflectance of the multispectral image were firstly taken as the input variables, and PLS, BPNN, and ELM algorithms were used to build the CNC prediction models, respectively.
The modeling results were listed in Table 6, while the scatter plots of the measured values and the predicted values of every model were also presented in Figure 6. As can be seen, the PLS model had the highest R v value of 0.73 and the lowest RMSE v value of 0.13 in both the calibration set and the validation set compared to the ELM and BPNN model. This indicated that PLS had better modeling performance. the first two VIs in Table 4, SRPI and NPCI, were recommended as efficient inputs. Furthermore, NGBDI, which had relatively higher grey relational degree with the CNC and a lower R with SRPI and NPCI, was also recommended as another efficient input variable.

Modeling with the Five-Band Reflectance
The five-band reflectance of the multispectral image were firstly taken as the input variables, and PLS, BPNN, and ELM algorithms were used to build the CNC prediction models, respectively.
The modeling results were listed in Table 6, while the scatter plots of the measured values and the predicted values of every model were also presented in Figure 6. As can be seen, the PLS model had the highest Rv value of 0.73 and the lowest RMSEv value of 0.13 in both the calibration set and the validation set compared to the ELM and BPNN model. This indicated that PLS had better modeling performance.

Modeling with VIs
In addition to the three recommended VIs (SRPI, NPCI and NGBDI), SIPI was also selected to build prediction models for comparison. Several models were established from a single VI or a combination of VIs with different modeling algorithms. The results were listed in Table 7. The PLS model still had a more accurate and balanced performance than the BPNN and ELM models. From the results of the single VI-based models, we could find that VIs which had higher GRD with the CNC had better performance in prediction models. From the results of the double and multiple VI-based models, the SRPI-and SIPI-based model showed lower accuracy than the SRPI-and NGBDI-based model, even though SIPI had a higher GRD than NGBDI. This proved the correctness of choosing NGBDI as the supplement variable rather than the others. The model with the highest accuracy (R v = 0.63) was established based on SRPI and NPCI and NGBDI, rather than all VIs, indicating that simply increasing the number of input variables does not necessarily improve the accuracy of the model. As long as the variables are selected correctly, fewer variables may bring higher accuracy to the model.

Modeling with the Five-Band Reflectance and VIs
Other types of models based on different combinations of the five-band reflectance and VIs were established, and the results were listed in Table 8. Of the different modeling algorithms, PLS still had the best performance among the three algorithms, BPNN had slightly lower accuracy than PLS, and ELM had an obvious lower accuracy than the other two algorithms. The PLS model based on FR and the three recommended VIs (SRPI, NPCI and NGBDI) had the highest accuracy, with the highest R v of 0.79 and the lowest RMSE v of 0.11. The scatter plots of the prediction result of the best model were presented in Figure 7. Higher accuracy based on the five-band reflectance combined with proper VIs was reached because this combination of input variables not only ensured the integrity of information, but also highlighted the spectral characteristic information. The results also demonstrated that the selection of input variables was crucial. 7. Higher accuracy based on the five-band reflectance combined with proper VIs was reached because this combination of input variables not only ensured the integrity of information, but also highlighted the spectral characteristic information. The results also demonstrated that the selection of input variables was crucial.

Irrigation Level Recognition
Sugarcane is a crop with high stem, large biomass, and long growth period. It has a large water requirement, as well as being dependent on fertilizers, especially in the early and middle growing stages (the seedling, tillering and elongating stage). Knowledge of water conditions during these stages is of great importance for sugarcane irrigation management.

Correlation Analysis
The correlations between the irrigation amounts and the spectral features, including the band reflectance and ten VIs, were firstly analyzed, as shown in Table 9. As can be seen, regardless of the amount of fertilizer applied, irrigation had high correlations of above 0.65 with red, blue, SRPI, NPCI and NGBDI. Red, blue and NPCI were negatively correlated with irrigation, while SRPI and NGBDI were positively correlated with irrigation. Among the five spectral bands, red, green and blue bands were negatively correlated with irrigation, while red edge and NIR bands were positively correlated with irrigation. This phenomenon is consistent with the spectral variation trend of the green plant (the healthier in growth condition, the lower the reflectance in visible range, and the higher the reflectance in red edge and NIR range).

Classification Results
Based on the five parameters of red, blue, SRPI, NPCI, and NGBDI, which had the highest correlations with irrigations SVM and BPNN were adopted to classify the irrigation levels. The results were listed in Table 10. As shown in Table 10, the OA of SVM for irrigation level recognition was 80.6%, while the OA of BPNN was only 61.7%. From the classification results of SVM, the UA and PA of irrigation_0 were the highest, reaching 94.9% and 87.0%, respectively. The accuracies of the other two levels generally exceeded 70%. This result reflects that the larger the difference in irrigation amount, the higher the classification accuracy. In addition, most of the misclassified samples were misclassified due to adjacent irrigation levels. For example, among the 108 samples in Irrigation_0, 94 samples were correctly identified, 13 samples were misidentified as Irrigation_180, and only one sample was misidentified as Irrigation_300. Among the total number of 63 misclassified samples, 62 samples were misclassified into adjacent levels, and only one sample was misclassified into a level far away. This indicated that the identification of irrigation levels using multispectral images had great potential and could be used for the recognition of crop water stress condition.

Discussion
The results in Tables 6-8 indicated that the PLS models had the best performance for sugarcane CNC prediction based on different input combinations, which is consistent with research into citrus CNC prediction by Liu et al. [63], grapevine LNC prediction by Moghimi et al. [64], etc. The sugarcane CNC prediction model showed a highest accuracy of R = 0.79 and RMSE = 0.11, which was also close to the previous studies by Liu et al. Compared to the five-band reflectance models in Table 6, VI-based models in Table 7 had obvious lower accuracy, indicating that taking VIs alone as the input would decrease the prediction accuracy comparing to taking the entire five-band reflectance as the inputs. It reflected that the VIs did not exert their characteristics, which should enhance the spectral features and reduce environmental interference [65]. The main reason was that VIs do have advantages when it comes to reducing the influence caused by uneven light (illumination) and different backgrounds. However, in this study, only one image of a small field was acquired in a very short time, meaning that the light difference and background difference was not significant. This made the contribution of VIs less than that of the whole spectral reflectance [66,67].
Regardless, VIs could still help to improve the modeling accuracy, and this was proved in the results listed in Table 8. Among different combinations of input spectral features, the five-band reflectance combined with the three VIs (SRPI, NPCI, and NGBDI) had the highest accuracy in CNC prediction, with the R v of 0.79 and the RMSE v of 0.11; this was 8.2% higher in R v , and 15.4% lower in RMSE v , than the five-band prediction model Moreover, all three of those VIs were calculated from the visible bands (SRPI and NPCI are both calculated from the green and red bands, while NGBDI is calculated from the blue and green bands), indicating that the visible bands contained more sensitive information for sugarcane CNC prediction. Ranjan et al. [68] explored the spectral characteristic of CNC in crops, with characteristic wavelengths mainly in the range of 430 nm, 460 nm, 640 nm, 910 nm, 1510 nm, 1940 nm, 2060 nm, 2180 nm, 2300 nm, and 2350 nm. In this study, the multispectral camera had a spectral range of about 450 nm-850 nm, which contained only the visible sensitive bands. This proves the rationality that the VI selected in this study is mainly concentrated in the visible range. As is generally known, different N inputs could lead to different leaf pigment concentrations, leaf internal structures, and canopy structures [64,69,70]. The visible bands are closely related to leaf pigments and canopy structures and, as such, offer a great potential for N prediction.
Furthermore, this research also discovered that the irrigation levels could be effectively classified based on the reflectance at red and blue, combined with the SRPI, NPCI and, NGBDI, which were the spectral features in the visible bands. Indeed, the NIR bands are generally more sensitive to plant water content. However, the most sensitive bands which sit between 1480 and 1500 nm are out of the range of the multispectral camera used in this study [69,70]. Insufficient water input could obviously affect plant metabolism, which indirectly affects the leaf pigment concentrations. Therefore, the irrigation level recognition model can achieve better classification accuracy by only using the three visible bands.
This study has achieved good results for CNC prediction and irrigation level classification based on multispectral remote sensing. Although research on crop monitoring based on UAV multispectral imagery have been widely carried out for more 10 years, but it is rarely applied in wide field management. Several bottlenecks need to be addressed at present. Sozzi et al. [4] compared the advantages and disadvantages of satellite-, plane-, and UAV-based multispectral imagery in variable rate N application in terms of cost, economic benefit, optical quality, and usage scenarios. It was pointed out that although satelliteand plane-based imagery have low optical quality and low resolution, they can provide applicable variable N rate suggestions, and bring economic benefits for large-scaled farms due to their relatively low cost. The UAV platform does have limits in acquisition cost and flight coverage at present. However, as the development of UAV technology and the increase requirement of UAVs, the cost can be significantly reduced and the battery performance can be enhanced in the future. With the emergence of automated UAV base stations and the reduction of image processing costs, the large-scale application of UAVs is just around the corner. By then, UAV remote sensing technology can be widely accepted in farm-scale crop monitoring with its flexible and autonomous acquisition style, high-quality image data, and low-cost validity.

Conclusions
High resolution multispectral images of a sugarcane field were collected by UAV, and its ability to predict CNC and irrigation levels was evaluated. The main conclusions were as follows: 1.
Ten VIs were used to determine sensitive spectral features for CNC and irrigation level prediction. The SRPI, NPCI, and NGBDI composed of the visible bands were sensitive to the sugarcane CNC as well as the irrigation level, and had a notable contribution to the accuracy improvement of the CNC prediction model.

2.
Different modeling algorithms based on different spectral features were compared in predicting sugarcane CNC. The PLS model had a clearly superior performance than the BPNN and ELM models. It was also crucial to select proper features among the band reflectance and VIs. The PLS model, based on the five-band reflectance combined with SRPI, NPCI, and NGBDI, had the highest R v of 0.79, and the lowest RMSE v of 0.11.

3.
Based on the correlation coefficients with irrigation levels, the red band, blue band, SRPI, NPCI, and NGBDI were adopted as the variables to classify irrigation levels based on SVM and BPNN, respectively. SVM reached an obviously superior performance compared to BPNN, with its overall classification accuracy of 80.6%.
High resolution multispectral images have been demonstrated as effective for CNC prediction and water irrigation level recognition. More adequate experiment could be conducted to collect more samples at different growth stages in sugarcane field. Time-or perioddependent prediction models could be studied to give users more accurate information.

Data Availability Statement:
The data presented in this study are available on request from the author (X.L.).