Estimation of Nitrogen in Rice Crops from UAV-Captured Images

Abstract

Leaf nitrogen (N) directly correlates to chlorophyll production, affecting crop growth and yield. Farmers use soil plant analysis development (SPAD) devices to calculate the amount of chlorophyll present in plants. However, monitoring large-scale crops using SPAD is prohibitively time-consuming and demanding. This paper presents an unmanned aerial vehicle (UAV) solution for estimating leaf N content in rice crops, from multispectral imagery. Our contribution is twofold: (i) a novel trajectory control strategy to reduce the angular wind-induced perturbations that affect image sampling accuracy during UAV flight, and (ii) machine learning models to estimate the canopy N via vegetation indices (VIs) obtained from the aerial imagery. This approach integrates an image processing algorithm using the GrabCut segmentation method with a guided filtering refinement process, to calculate the VIs according to the plots of interest. Three machine learning methods based on multivariable linear regressions (MLR), support vector machines (SVM), and neural networks (NN), were applied and compared through the entire phonological cycle of the crop: vegetative (V), reproductive (R), and ripening (Ri). Correlations were obtained by comparing our methods against an assembled ground-truth of SPAD measurements. The higher N correlations were achieved with NN: $0.98$ (V), $0.94$ (R), and $0.89$ (Ri). We claim that the proposed UAV stabilization control algorithm significantly improves on the N-to-SPAD correlations by minimizing wind perturbations in real-time and reducing the need for offline image corrections.

1. Introduction

Farmers use soil plant analysis development (SPAD) devices as a field diagnostic tool to estimate nitrogen (N) content in the plant and to predict grain yield [1,2]. For rice crops, N provides critical insight into plant-growth, crop yield, and biomass accumulation [3]. Furthermore, monitoring N allows farmers to properly adapt crop management practices [4,5]. However, using SPAD devices for crop diagnosis is still time consuming. With the advent of low-cost unmanned aerial vehicles (UAVs), several authors have reported faster and accurate remote sensing tools and methods [6,7] to estimate the N canopy from aerial multispectral imagery [8,9]. One of the most common methods used, relies on sensing the canopy reflectance in both visible (VIS) and near-infrared (NIR) wavelengths, using hyperspectral sensors [10,11,12,13].

Multiple vegetation indices (VI) have been established to associate specific spectral reflectance wavelengths with the crop variables of interest [14,15]. Although these VIs can be used to estimate the N contents in rice plants, there is no single standard set of VIs that works across all crop stages, rice varieties, soils, and atmospheric conditions. Most of the existing body of work associated with the estimation of N contents in plants, has been conducted using multi-variable regressions [5,16], in which linear functions are defined by combining and weighting each VI according to the regression model. This approach is simple, but sometimes inaccurate since the evolution of the crop tends to exhibit highly nonlinear effects that affect the N content over time. Machine and deep learning methods have recently gained traction, in solving these issues. Machine learning has the potential to drive significant breakthroughs in plant phenotyping [17,18,19,20].

For canopy N estimation, training machine learning algorithms requires the precise extraction of VI features from the aerial multispectral imagery. Several authors have applied data fusion methods from different sensors [21,22,23] for applying image mosaicing techniques [24,25,26,27], computing crop surface models based on dense image matching (DIM) techniques [28], or applying photogrammetry methods that are computationally expensive [26,27,29,30]. Other approaches rely on the real-time segmentation and image registration for the extraction of relevant features associated with the leaf/canopy N, including edge detection thresholding, color histograms, and clustering (otsu, K-means, watershed) [31,32,33].

In this work, we improve on the stabilization control of the UAV for acquiring accurate plot imagery, instead of relying on mosaicing or photogrammetry methods. Commercial UAV quadrotors used for the monitoring of the canopy N usually come with a proportional-integral-derivative (PID) waypoint navigation autopilot. The lack of proper UAV angular stabilization in the presence of large wind perturbations limits the accuracy of image registration algorithms, therefore affecting the estimation of canopy N content through image processing. To address this problem, we demonstrate a nonlinear trajectory-tracking controller that enables precise UAV positioning through a stabilization/attitude control loop. We call this method backstepping+ desired angular acceleration function (BS + DAF). This approach incorporates UAV aerodynamics information within the control law to produce roll and pitch acceleration commands to reduce abrupt angular acceleration changes caused by external wind disturbances. This added compliance enables the UAV to hover over the crop plots precisely, which in turn allows for capturing imagery that can be individually processed in real-time.

Here, we combine an autonomous UAV to acquire multispectral imagery from a crop, with machine learning methods as a means to high-throughput nitrogen content estimation. Three major contributions are involved:

(i)

The development and validation of a novel UAV attitude control called BS+DAF. The UAV captures multispectral imagery with relevant above-ground data that must correlate with the SPAD measurements at the ground-level of the crop. The proposed controller is aimed at maintaining precise angular stabilization during flight by properly rejecting wind disturbances, enabling improvements in the estimations of N.

(ii)

The application and validation of a segmentation algorithm called GrabCut [34], instead of the traditional edge detection thresholding methods, color histograms, and clustering techniques used in agriculture. The GrabCut approach achieves smooth pixel information with richer detail of the canopy structure, enabling the accurate segmentation of the multispectral imagery and the proper VI-based feature extraction.

(iii)

The integration of machine learning methods to process nonlinear N dynamics with the calculations of the VIs during all stages of crop growth. A comprehensive characterization of the crop by designing a ground-truth dataset with several contrasting rice genotypes and accurate direct measurements of leaf nitrogen (training model).

2. Materials and Methods

2.1. System

Figure 1 details the proposed approach. Our UAV is a quadrotor called the AscTec Pelican (manufactured by Intel’s Ascending Technologies (https://robots.ros.org/astec-pelican-and-hummingbird/)). This UAV comes with a C++ software development kit (SDK) that enables onboard code development with ease. A high-level Atom Intel processor (HLP) offers enough computing power to run solutions in real-time, whereas a low-level processor (LLP) handles the sensor data fusion and rotor control with an update rate of 1kHz. As shown in Figure 1, we developed a closed-loop attitude controller to drive the low-level rotor’s control running in the LLP. This control method depends on the dynamics of the UAV to properly reject wind disturbances and keep the UAV steady. During flight, a dataset of multispectral imagery is precisely collected aiming at the above-ground estimation of nitrogen by using machine learning methods. In previous work [35], our UAV system was tested during two years of in-field experiments with the aim of estimating above ground biomass accumulation. Here, we have extended the capabilities of the system in Figure 1 by developing and integrating novel methods for UAV flight control, multispectral imagery segmentation for VI feature extraction, and their impact on nitrogen estimation using machine learning algorithms.

The UAV was equipped with the Parrot Sequoia multispectral sensor (https://www.parrot.com/business-solutions-us/parrot-professional/parrot-sequoia) fabricated with 4 separate bands: green, red, red-edge, and near-infrared. The camera offers a resolution of $1280\times 960$ for each independent spectral sensor, yielding a crop-to-image resolution of 1.33 cm/pixel according to the flying altitude of 20 m. In addition, the multispectral camera comes with a radiometric calibration target that enables offline reflectance calibration across the spectrum of light captured by the camera, and an onboard irradiance sensor designed to correct images for illumination differences, all of which enables its outstanding performance in cloudy conditions, as evaluated in [36]. In our application, this calibration procedure was executed prior to any flight mission of the UAV. Additional sensors such as an IMU, magnetometer, and GPS are also embedded within the sunshine sensor. Figure 2a depicts the mounted camera setup, while

Table 1. General specifications.

Specification	Description
UAV
UAV Processing	High-level: Intel®CoreTM i7|Low-level: Intel Atom 1.6 GHz
Flight altitude over ground	20 m
Flight linear speed	1.5 ms${}^{-1}$
Autonomy	25 min
NAVIGATION
DGPS Piksi GNSS Module ${}^1$	HW version: 00108-10
DGPS Processing	dual-core ARM Cortex-A9 at 666 MHz
Real-time kinematic positioning (RTK)	10 Hz: GPS, GLONASS, BeiDou, Galileo
DGPS Position accuracy	35 cm
IMU rate and accuracy	1 KHz|Heading: 0.2${}^{\circ }$, Pitch/Roll: 0.03${}^{\circ }$
IMAGERY
Parrot Sequoia multispectral camera	Spectral bands resolution: $1280\times 960$ pixels. Green (530–570 nm), Red (640–680 nm), Red-Edge (730–740 nm), Near-Infrared (770–810 nm). RGB sensor resolution: $4608\times 3456$ pixels.
Imagery resolution at 20 m altitude (Ground Sampling Distance—GSD) ${}^2$	RGB: 0.37 cm/pixel. Spectral-bands: 1.33 cm/pixel
Onboard irradiance sensor	real-time illumination image correction
CROPS
Experiments location	Lat: 4${}^{\circ }{1}^{\prime }{37.85}^{″}$ N, Lon: 73${}^{\circ }{28}^{\prime }{28.65}^{″}$ W International Center for Tropical Agriculture (CIAT), Department of Meta, Colombia ${}^3$
Weather conditions ${}^4$	Dry season (June-September): Average temperatures: L: 22${}^{\circ }$–H: 31 ${}^{\circ }$C with average cloudiness of 30%.

${}^1$ Differential Global Positioning System (DGPS) Piksi module: https://www.swiftnav.com/piksi-multi. ${}^2$ Ground Sampling Distance (GSD) calculator: https://support.pix4d.com/hc/en-us/articles/202560249-TOOLS-GSD-calculator#gsctab=0. ${}^3$ International Center for Tropical Agriculture (CIAT) website: https://ciat.cgiar.org. ${}^4$ Weather information: https://weatherspark.com/y/24273/Average-Weather-in-Villavicencio-Colombia-Year-Round.

details the general specifications of our system.

2.2. Rice Crops

The crops were designed with 3 spatial repetitions containing 8 contrasting rice genotypes in terms of N concentration, biomass accumulation, and flowering: FED50, MG2, ELWEE, NORUNKAN, IR64, AZUCENA, UPLRI7, and LINE 23. These rice varieties have a phenological cycle ranging between 86–101 days , as detailed in Figure 3b. The experimental design consisted of six months of in-field sampling with three different planting crops. For each experiment, the amount of applied nitrogen and water was modified as follows:

• Experiment 1: 200 kgha${}^{-1}$ of nitrogen and limited crop irrigation.
• Experiment 2: 200 kgha${}^{-1}$ of nitrogen with permanent crop irrigation.
• Experiment 3: 100 kgha${}^{-1}$ of nitrogen with permanent crop irrigation.

Experimental trials were conducted during the dry season. To assess the effects of limited and permanent irrigation on the crops, leaf temperature (MultispeQ (https://www.photosynq.com)), and soil humidity (AquaPro (http://aquaprosensors.com)) were constantly monitored. Fertilizers were applied accordingly in order to maintain the crops through the phenological cycle. Given that, the same amount of fertilizers were applied for the experiments: 60 kgha${}^{-1}$ of $P_2{O}_5$ (diammonium phosphate) and 130 kgha${}^{-1}$ of $K_{2}O$ (potassium chloride), as detailed in the experimental protocol available in the Supplementary Materials section.

Figure 2 details some characteristics of the crop. As shown in plots (a), (b), the length of a plot was about 2.75 m with a distance between plants of 25 cm and 30 cm between rows. Within each plot, we defined 6 sampled areas with 1 linear meter in length (containing four plants). A ground-truth was defined based on the direct measurements of plant chlorophyll using a SPAD 502 Plus meter (Konica-Minolta) over these sampled areas. Measurements from the crop were obtained during three stages of rice growth: vegetative, reproductive, and ripening. Figure 3b details the specifications of the crop experiments reported here.

Table 2. An example of a ground-truth dataset. The crop field was designed with three spatial repetitions (Rep) containing 8 contrasting rice genotypes.

Plot	Genotype	Rep	SPAD
1	AZUCENA	1	56.55
2	ELWEE	1	47.80
3	LÍNEA 23	1	54.55
4	UPLRI7	1	46.32
5	NORUNKAN	1	43.30
6	IR64	1	32.91
7	FED50	1	47.06
8	MG2	1	43.36
9	AZUCENA	2	49.26
10	IR64	2	42.59
11	LÍNEA 23	2	49.82
12	UPLRI7	2	48.15
13	ELWEE	2	41.29
14	FED50	2	46.29
15	NORUNKAN	2	42.82
16	MG2	2	38.96
17	FED50	3	49.20
18	UPLRI7	3	40.88
19	IR64	3	40.23
20	AZUCENA	3	55.81
21	NORUNKAN	3	43.96
22	ELWEE	3	49.70
23	LÍNEA 23	3	46.45
24	MG2	3	42.67

details how the dataset was collected, in which the SPAD value corresponds to the average of 24 measurements conducted over each plot, as depicted in Figure 3b. All ground-truth data are available through the Open Science Framework in the Supplementary Materials section. On the other hand, Figure 3a shows the linear correlation used to relate the measured SPAD value with the leaf-blade N concentration [37].

As mentioned, vegetation indices are widely used to quantify both plant and soil variables by associating certain spectral reflectances that are highly related to variations in canopy chemical components such as nitrogen. From the extensive list of vegetation indices (VIs) available [38,39,40,41], we selected 7 VIs with sufficient experimental evidence and quantitative trait loci (QTL)-based characterization regarding their high correlation with rice nitrogen [42].

Table 3. Selected near infrared vegetation indices (extracted features). The ${\rho }_f$ term denotes the reflectance of the frequency f).

Name	Equation
Normalized Difference Vegetation Index [38]	$\mathrm{NDVI}=\frac{{\rho }_{780}-{\rho }_{670}}{{\rho }_{780}+{\rho }_{670}}$
Green Normalized Difference Vegetation Index [39]	$\mathrm{GNDVI}=\frac{{\rho }_{780}-{\rho }_{500}}{{\rho }_{780}+{\rho }_{500}}$
Simple Ratio [38]	$\mathrm{SR}=\frac{{\rho }_{780}}{{\rho }_{670}}$
Soil-Adjusted Vegetation Index [39,40]	$\mathrm{SAVI}=(1+\mathrm{L})\left(\frac{{\rho }_{800}-{\rho }_{670}}{{\rho }_{800}+{\rho }_{670}+\mathrm{L}}\right)$ with $\mathrm{L}=0.5$
Modified SAVI [40]	$\mathrm{MSAVI}=\frac{1}{2}\left(2{\rho }_{800}+1-\sqrt{{(2{\rho }_{800}+1)}^2-8({\rho }_{800}-{\rho }_{670})}\right)$
Triangular Vegetation Index [39]	$\mathrm{TVI}=\frac{1}{2}\left(120({\rho }_{780}-{\rho }_{500})-200({\rho }_{670}-{\rho }_{500})\right)$
Corrected Transformed Vegetation Index [42]	$\mathrm{CTVI}=\frac{\mathrm{NDVI}+0.5}{|\mathrm{NDVI}+0.5|}\sqrt{|\mathrm{NDVI}+0.5|}$

details the selected VIs. These formulas are applied in different wavelengths for taking into account the changes in canopy color, since several factors affect the spectral reflectances of the crop: solar radiation, plant morphology and color, leaf angles, undergrowth, soil characteristics, and water.

In our system, the parrot sequoia camera has a solar radiation sensor that compensates the light variations in the canopy. The change in the rice canopy color is perhaps the most notable variation during the phenological cycle. In the vegetative stage, the green color is predominant whereas in the reproductive stage, panicle formation starts and thus yellow features appear in the images. In ripening, the maturation of the plants occur while the leaves begin to senesce, turning the yellow color predominant. These changes can be observed in Figure 2c. Furthermore, different wavelengths of light have a different level of plant absorption depending on the leaf composition given by several genetic traits. In particular, the relation between the selected VIs in Table 3 with the photosynthetic activity and canopy structural properties has allowed the association of certain spectral reflectances that are highly related to the physico-chemical canopy N variations in plants, especially the green, red, red-edge, and near infrared bands.

The selected VIs exhibit a strong dependence on the NIR reflectance due to leaf chlorophyll absorption, providing an accurate approach to determine the health status of the plants and the canopy N. Most of the existing body of research focused on multispectral-based N estimations [14,21,40,42,43], combine the information provided by several vegetation indices in order to avoid saturation issues. For instance, the NDVI, which is one of the most common VIs used, tends to saturate with dense vegetation. In turn, the NDVI alone is not accurate during the reproductive and ripening stages of rice growth. Here, we combine several VIs across the crop stages, to ensure data on wavelengths located in the red-edge and another spectral reflectances that accurately express the healthy status of the leaves (higher NIR and green-band readings).

2.3. UAV Stabilization Control

The estimation of the canopy N requires the precise extraction of Vegetative Index features from the acquired multispectral imagery. Capturing aerial images during UAV hovering that accurately matches with the GPS positions of the SPAD measurements at ground-level is essential.

As previously detailed in Figure 1, three closed-loop controllers are needed to regulate: (i) the X-Y position based on GPS feedback, (ii) the Z altitude based on barometric pressure and laser readings (pointing downwards), and (iii) the $\varphi ,\theta ,\psi$ attitude based on IMU data. The UAV is constantly subjected to wind disturbances that cause unsteady angular motions and therefore imprecise trajectory tracking. Consequently, aerial imagery captured across the crop with the UAV system is affected. To overcome this issue, we replaced both roll and pitch PID-based controllers by a robust nonlinear backstepping (BS) control.

The classical BS method has several advantages. It explicitly takes into account the nonlinearities of the UAV model defined in Equation (A2) (see Appendix A), and most importantly, allows the incorporation of a virtual control law to regulate angular accelerations. For this, we have derived a desired acceleration function (DAF) for roll and pitch. This enhanced controller is called backstepping+DAF (BS+DAF). Our goal is to use the dynamics equations of motion (EoM) defined in Algorithm A1 within the control law in order to compensate for abrupt angular acceleration changes, concretely in roll and pitch. The DAF terms improve the responsiveness of the control law to external perturbative forces. Appendix B details the control law derivation for roll and pitch.

The BS + DAF control supports on the Lyapunov stability concept that guarantee asymptotic stabilization around the equilibrium points. For our application, we require both roll $\varphi$ a pitch $\theta$ angles to remain in zero while the UAV is hovering above the crop for capturing multispectral imagery, i.e., $e_{\varphi }={\varphi }^d-\varphi \to 0$ and $e_{\theta }={\theta }^d-\theta \to 0$. Otherwise, the set-point references for both roll and pitch controllers are defined by the X-Y position controller.

As previously mentioned, high wind-speed perturbations affect the UAV during hovering. Nonetheless, our controller law is sensitive to angular acceleration changes caused by external wind disturbances, thanks to both error dynamics ($e_2$) and the DAF terms (${\ddot{\varphi }}^d$) in Equation (A18) (Appendix B). Given that, we present how to model the effects of the wind over the UAV aerodynamics, by following the wind effect model (WEM) developed by [44]. Figure 4 details the control architecture and how the WEM module has been incorporated as an external disturbance acting on the UAV body frame.

The WEM defines the relationship between the oblique airflow acting on a propeller and the resulting forces and moments applied onto the UAV. In our quadrotor system, the oblique airflow form and angle of 90^o with the propeller azimuth angle. In terms of aerodynamics, the incoming airflow induces an increase in vertical thrust ($T_{oi}$), a side force (F${}_{wind}$), and a pitching moment (M${}_{wind}$). The WEM presented by [44] uses blade element momentum theory to determine the aforementioned effects. Here, we are interested in incorporating the effects caused to the thrust and the pitching moment (M${}_{wind}$) as a function of the magnitude of a given wind speed vector $V_{wind}$ and the rotational speed $\omega$ of each propeller driven by the rotor’s model. In this regard, the WEM is defined as:

$$\begin{array}{c}{M}_{wind}=\left[\begin{array}{ccc}cos\left(\beta +\frac{\pi }{2}\right)& sin\left(\beta +\frac{\pi }{2}\right)& 0\\ -sin\left(\beta +\frac{\pi }{2}\right)& cos\left(\beta +\frac{\pi }{2}\right)& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}{M}_{prop}\\ 0\\ 0\end{array}\right]+\frac{l}{2}\left[\begin{array}{c}{C}_4{T}_4-C_2{T}_2\\ C_1{T}_1-C_3{T}_3\\ 0\end{array}\right]\hfill \end{array}$$

(1)

The angle $\beta$ determines the direction of the $V_{wind}$ vector. Depending on the direction of the applied relative wind, certain propellers will be more affected by this external disturbance. Tran et al., in [44,45], propose the use of a weight function that assigns values between 0 to 1, where the maximum value of 1 means the propeller is directly exposed to the wind. In Equation (1), C corresponds to the weighting vector that affects the thrust T generated by each propeller. The parameter l is the distance between non-adjacent propellers and $M_{prop}$ is:

$$\begin{array}{c}{M}_{prop}=\left[\begin{array}{cccc}{C}_1& C_2& C_3& C_4\end{array}\right]\left[\begin{array}{c}{\tau }_1\\ {\tau }_2\\ {\tau }_3\\ {\tau }_4\end{array}\right]\hfill \end{array}$$

(2)

Therefore, the scale of the weighting vector C is defined based on the magnitude of the applied wind speed vector and the corresponding direction (the $\beta$ angle between the relative wind vector pointing towards the UAV body frame). In addition, the rotor’s speeds $\omega$ are required to calculate the augmented thrusts T and torques $\tau$ (weighted by C). Further details on this model can be found in [44,45].

Section 3.1 and Section 3.3 present the results regarding the impact of precise UAV positioning on the accurate estimation of the canopy N through the entire phenological cycle.

2.4. Multispectral Image Segmentation

Multispectral imagery is evaluated by using a multispectral image segmentation method called GrabCut [34]. The original method is widely used for its ease of implementation and for the excellent results in generating a binary classification; however, it suffers from the drawback of being a semi-manual algorithm. The original GrabCut method requires a manual input during the algorithm iteration in order to properly determine both background and foreground pixel values.

This section introduces a modified version of the GrabCut algorithm that is fully automatic, thanks to an added refinement step using a guided filtering [46] to extract the relevant pixel information associated with the plant’s canopy. Our approach solves an optimization problem using an energy function that allows the proper labeling of texture in the multispectral image by using a Gaussian mixture model. As mentioned, we use a refinement process based on a guided-filter by taking into account information from all band channels of the multispectral camera: green, red, red-edge, and near-infrared. The resultant multispectral image-mask includes only relevant pixel information that accurately represents the canopy for the estimation of N.

2.4.1. GrabCut Segmentation

The GrabCut method requires an initial model known as a trimap. This model consists of a partition of the input image into three regions: a foreground $T_F$, a background $T_B$, and a region with pixels that result from the combination of the foreground and background colors $T_U$.

The image is modeled as an array of values $\mathbf{z}=(z_1,\dots ,z_N)$, and the output segmentation can be a channel of values between 0 and 1 or a hard segmentation with a binary assignment (0 or 1). The segmentation is written as $\underline{\alpha }=({\alpha }_1,\dots {\alpha }_N)$ with $0\le {\alpha }_n\le 1$, or ${\alpha }_n=\{1,0\}$.

The GrabCut algorithm also requires a model for the distribution of the foreground/background colors and gray levels. This distribution can be represented as a histogram directly assembled from $T_F$ and $T_B$, as: $\underline{\theta }=\{h(z;a),a=0,1\}$. Under this trimap model, the segmentation algorithm determines the values of $\alpha$ from the image $\mathbf{z}$ and the distribution model $\underline{\theta }$. The $\alpha$ values are calculated from an energy minimization function, as:

$$\mathbf{E}(\underline{\alpha },\underline{\theta },\mathbf{z})=U(\underline{\alpha },\underline{\theta },\mathbf{z})+V(\underline{\alpha },\mathbf{z}),$$

(3)

where the sub-function U evaluates the fitness by assigning a low score if the segmentation ($\alpha$) is consistent with the image $\mathbf{z}$, defined as follows:

$$U(\alpha ,\theta ,\mathbf{z})=\sum _n-lnh(z_n;{\alpha }_n)$$

(4)

Instead of using the histogram as the estimator of the probability density function, the algorithm uses a Gaussian mixture model in order to take into account the information coming from all channels.

$$V(\underline{\alpha },\mathbf{z})=\gamma \sum _{(m,n)\in \mathbf{C}}[{\alpha }_n\ne {\alpha }_m]e^{-\beta {(z_m-z_n)}^2}$$

(5)

In Equation (5), the sum set C $\in 3\times 3$ refers to the neighbors pixels in a given window, and the the term $\beta$ can be calculated according to Equation (6).

$$\beta =\frac{1}{\left(2E\left|\langle {(z_m-z_n)}^2\right|\rangle \right)}$$

(6)

By using the global minimum in 7, the image segmentation is estimated as:

$$\underline{\widehat{\alpha }}=arg\underset{\underline{\alpha }}{min}\mathbf{E}(\underline{\alpha },\underline{\theta },\mathbf{z})$$

(7)

Algorithm 1 details the original GrabCut method. In order to eliminate the third manual step of the algorithm, a fixed background image $T_B$ mask is used during the iteration, and a guided-filter refinement process is applied to achieve an automatic segmentation, as detailed in Algorithm 2.  

Algorithm 1: Original GrabCut algorithm.

Step 1: Initialization of $T_B$, $T_F$, $T_U$. Initialize the trimap $T_F$ as the foreground Initialize the trimap $T_B$ as the background All remaining pixels are set as a possible foreground pixels $T_{UF}$ Step 2: Iterative minimization. Assign and learn GMM Minimize energy function Estimate segmentation image $\alpha$ IF image $\alpha$ have visible errors GO to Step 3 Step 3: User editing. Use the segmented image $\alpha$ as the new possible foreground pixels $T_{UF}$ Input pixel hints for $T_B$ and $T_F$ GO to Step 2

Algorithm 2: Modified GrabCut with guided filter calculation. In the algorithm, $E_r\left(I\right)$ denotes a function that calculates the image mean over a radius r, the operations $.*$ and $./$ denotes the matrix element-wise calculation, and q is the image output.

Step 1: Initialization of $T_B$, $T_F$ and $T_U$. Initialize the trimap $T_F$ as the foreground Initialize the trimap $T_B$ as the background All remaining pixels are set as a possible foreground pixels $T_{UF}$ Step 2: Iterative minimization. Assign and learn GMM Minimize energy function Estimate segmentation image $\alpha$ Step 3:Input image p, input guidance I, radius r, and regularization $ϵ$. 1: ${\mu }_I\leftarrow E_r\left(I\right)$, ${\mu }_p\leftarrow E_r\left(I\right)$, $Cor{r}_I\leftarrow E_r(I.*I)$, $Cor{r}_{Ip}\leftarrow E_r(I.*p)$. 2: ${\sigma }_I^2\leftarrow Cor{r}_I-{\mu }_I$.*${\mu }_I$, ${\sigma }_{Ip}^2\leftarrow Cor{r}_{Ip}-{\mu }_I$.*${\mu }_p$ 3: $a\leftarrow {\sigma }_{Ip}^2$./$({\sigma }_I^2+ϵ)$, $b\leftarrow {\mu }_p-a.*{\mu }_I$ 4: ${\mu }_a\leftarrow E_r\left(a\right)$, ${\mu }_b\leftarrow E_r\left(b\right)$ 5: $q={\mu }_a$.*$I+{\mu }_b$

2.4.2. Guided Filter Refinement

The guided filter (GF) [46] can be defined as a convolutional Bilateral Filter with a faster response due to its $O\left(N\right)$ computational complexity. In this regard, the output of each pixel denoted as q can be expressed as a weighted average over the convolutional window W ($i,j$ denote the pixel coordinates):

$$q_i=\sum _j{W}_{ij}{p}_j$$

(8)

The GF implies that an image can be filtered using the radiance of another image as guidance. We exploited this concept to filter our segmented plot mask created with the GrabCut algorithm, with the aim of refining the segmented image with the original image as guidance. In this regard, the weight used by the GF is given by:

$$W_{ij}^{GF}\left(I\right)=\frac{1}{{\left|\omega \right|}^2}\sum _{k;(i,j)\in {\omega }_k}\left(1+\frac{(I_i-{\mu }_k)(I_j-{\mu }_k)}{{\sigma }_k^2+ϵ}\right),$$

(9)

where $W_{ij}^{GF}$ depends entirely of the guidance image I. The parameters ${\mu }_k$ and ${\sigma }_k^2$ are the mean and variance of the guidance image I estimated over a window $w_k$, $ϵ$ denotes a regularization parameter and $\left|\omega \right|$ is the number of pixels in the window $w_k$.

Section 3.2 presents the results of applying the proposed modified version of GrabCut, by comparing the method against traditional segmentation techniques such as thresholding, K-means, meanshift, but also against the original semi-manual GrabCut method.

2.5. Machine Learning for N Estimations

As detailed in Table 2, the ground-truth for training machine learning (ML) algorithms was defined based on the direct measurements of plant chlorophyll using a SPAD 502 Plus meter (Konica-Minolta) over these sampled areas, as depicted in Figure 1. Datasets contain the measured SPAD value that directly correlates with the leaf-blade N concentrations by following the linear correlation [37] defined in Figure 3a. In this regard, measurements from the crop were obtained during three stages of rice growth: vegetative, reproductive, and ripening, in which 3 trials were conducted per crop stage. These datasets are the result of 10 flights per trial, capturing around 500 images, and yielding a database of $1500$ images per stage. Since 3 trials were conducted per crop stage, around $13,500$ images were processed in this work. Figure 3b details the experimental specifications.

The ML methods were trained with a set of images accounting for the $60\%$ of the entire database. For the final estimations of leaf nitrogen, we used the remaining $40\%$ of the database (testing phase of the ML methods). The entire imagery dataset and the ground-truth available in the Supplementary Materials section. Here, we report on the use of classical ML methods based on multi-variable linear regressions (MLR), support vector machines (SVM), and artificial neural networks (NN) for the estimation of the canopy N.

For MLR models, the VIs from Table 3 were combined by following the formula: ${\displaystyle N={\beta }_0+\sum _{k=1}^7{\beta }_k{\left(VI\right)}_k}$, where the parameter ${\beta }_k$ changes accordingly to the growth stage of the plant, by weighting each VI independently. In this regard, each crop stage will have an independent MLR model that linearly fits the N content.

With the aim of improving the adaptability of the estimation models by considering several nonlinear effects (dense vegetation, optical properties of the soil, and changes in canopy color), this section presents the estimation results using support vector machines (SVM) and artificial neural networks (NN). SVM models were used with different kernel functions with the aim of determining the proper mathematical function for mapping the input data. Six different kernels based on linear, quadratic, cubic, and Gaussian models will be tested for each crop stage independently.

Contrarily to the MLR and SVM methods, in which N estimation models are determined for each crop stage independently, neural networks enable the combination of the entire dataset (SPAD values and VI extracted features) into a single model that fits for the entire crop phenological stages. Several non-linear training functions are tested with different hidden layers. In addition, we discarded the use of deep-learning methods such as convolutional neural networks (CNN), due to the high computational costs associated to the pooling through lots of hidden layers in order to detect data features. For this application, we use well-known vegetative index features (reported in Table 3) that have been widely used and validated in the specialized literature [14,40,42]. Other image-based features such as color, structure, and morphology do not work well with multispectral imagery of 1.2 mega-pixel in resolution, compared to the 16 mega-pixel in the RGB image. In fact, the main advantage of using VIs (as features for training), relies on having information at different wavelengths, providing key information of the plant health status related to N.

In Section 3.3, we report a comprehensive comparison among multi-variable linear regressions (MLR), support vector machines (SVM), and artificial neural networks (NN) for the estimation of the canopy N. We used three metrics based on the root mean square error ($RMSE$), Pearson’s linear correlation (r), and coefficient of determination ($R^2$), detailed as follows:

$$\begin{array}{c}RMSE=\sqrt{\frac{1}{n}{\sum }_{i=1}^n{({\widehat{N}}_i-N_i)}^2}\hfill \\ \\ r=\frac{{\sum }_{i=1}^n(N_i-\overline{N})}{\sqrt{{\sum }_{i=1}^n{(N_i-\overline{N})}^2}}\hfill \\ \\ R^2=\frac{{\sum }_{i=1}^n{(N_i-{\widehat{N}}_i)}^2}{{\sum }_{i=1}^n{(N_i-\overline{N})}^2}\hfill \end{array}$$

(10)

where n is the total number of samples, N is the measured value (SPAD scale cf. Figure 3a), $\widehat{N}$ is the estimated nitrogen, and $\overline{N}$ is the mean of the N values. In this application, the Pearson’s metric indicates the linear relationship between the estimated value of N VS the measured one (ground-truth). On the other hand, the $R^2$ metric is useful since it penalizes the variance of the estimated value from the measured one.

As mentioned, the ML models require the calculation of the VI formulas presented in Table 3, in order to determine the input feature vector. Every image captured by the UAV is geo-referenced with the DGPS system with a position accuracy of 35 [cm] (see Table 1). Given that, aerial imagery is registered by matching consecutive frames according to a homography computed with the oriented features from accelerated segment test (FAST) and rotated binary robust independent elementary feature (BRIEF) computer vision techniques [47]. The UAV path is planned by ensuring an image overlapping of $60\%$, which allows for the precise matching with the GPS coordinates of the ground-level markers to ensure a proper comparison between the aerial-estimations and ground-measurements of plant N. Part of this procedure is described in previous work reported in [48]. In turn, the metrics described in Equation (10) report on the performance of the ML models based on the aerial-ground data matching of canopy N. The aforementioned geo-referencing process is conducted by using the affine transformation, a 1st order polynomial function that relates the GPS latitude and longitude values with the pixel coordinates within the image. This procedure is also detailed in previous work reported in [48].

Figure 5 details the procedure required by the machine-learning models. After registration, images are segmented by using the proposed GrabCut method. Although all pixels in the image are evaluated in Algorithm 2, we use pixel clusters of $10\times 10$ to calculate the VI formulas, as shown in Figure 5a. After segmentation, pixels representing the rice canopy are separated from the background, as shown in Figure 5b. All vegetation indices are calculated as the average within each image sub-region.

At canopy-level, several factors affect the spectral reflectances of the crop: solar radiation and weather conditions, plant morphology and color, leaf angles, undergrowth, soil characteristics, and ground water layers. As mentioned in Section 2.1, the multispectral camera comes with an integrated sunshine sensor to compensate light variations in the resultant image. In addition, the GrabCut segmentation method deals with the filtering of undergrowth and other soil noises. Despite that, it remains crucial to validate the accuracy of the selected VIs, since the estimation of N depends on the accuracy and reliability of the extracted features. In this regard, Figure 5c presents several VI calculations in order to analyze the variance of the VI features through the entire phenological cycle of the crop. For this test, we calculated the VI formulas from 360${}^{\circ }$ random images per crop stage under different environmental conditions. We show the most representative VIs to nitrogen variations: simple ratio (SR), normalized difference vegetation index (NDVI), green normalized difference vegetation index (GNDVI), and soil-adjusted vegetation index (MSAVI). As observed, the low variance exhibited through the entire phenological cycle allows us to use the VI features as inputs to the ML models detailed in Figure 5d.

3. Results

The experiments reported in this paper were carried out during 2018 in the rice farms of the Center of International Agriculture (CIAT). Trials were conducted during the dry season (June-September), as detailed in Table 1. We conducted 3 trials per crop stage (vegetative, reproductive, and ripening).

3.1. UAV Control Results

Figure 6 shows simulation results to evaluate the performance of the proposed controller in terms of wind-disturbance rejection. In plot (a), the desired trajectory (red line) was defined for the UAV to cover the crop while following the trapezoidal velocity profile shown in plot (b). This trajectory profile enables the UAV to hover at certain points (black dots) to capture aerial images. The maximum UAV velocity was set to 1.5 ms${}^{-1}$.

In plot (a), a wind disturbance ($V_{wing}=9$ ms${}^{-1}$) was added at the starting point of the trajectory, causing the UAV to mismatch the desired path (the $V_{wing}$ vector was applied onto the $y+$ direction). The response of the UAV driven by the BS-DAF attitude controller corresponds to the black line, whereas the PID controller is the green line. As observed, the BS + DAF immediately counteracted the disturbance by generating the corresponding roll command $\varphi$ depicted in plot (d). This response allowed the UAV to follow the path precisely, maintaining the position error along the y axis near to zero, as shown in plot (c). Under this scenario, our system obtained a maximum tracking error of $2\%$.

A second wind disturbance was added when the UAV reached 10 m in altitude. In this case, $V_{wing}$ was applied onto the z-direction (z points downwards). As observed, an augmented thrust caused a large altitude tracking error that affected the UAV for both control schemes similarly. Unlike for roll and pitch, the BS+DAF control does not drive the altitude loop. In general, the UAV was able to position in the hovering knot-points accurately (black dots in plot a), maintaining a minimum error with the GPS waypoints. For the rest of the trajectory, the results demonstrate that the backstepping + DAF is accurate and reliable to reject external wind disturbances.

3.2. Image Segmentation Results

Algorithms 1 and 2 are applied for each input image in an iterative manner. In Figure 7a, two markers placed as geo-referenced guiding points appear in the RGB image (red dots in the top-right corner). The algorithm provides to the user with the option to manually select these points in order to remove them from the segmentation process. An example of semi-manual GrabCut segmentation with minimal user interaction is shown in Figure 7d. If those points are not manually removed, the corresponding pixels will be classified either in the background of foreground cluster automatically. Subsequently, plots (b) and (c) show how the initial foreground and background are defined. Plot (d) shows an example of applying Algorithm 1 to the input image from (a), whereas plot (e) shows the results when combining the GrabCut with the guided-filter (GF) approach, yielding an automatic segmentation process, as detailed in Algorithm 2. As shown in plot (e), the GF refinement achieves richer detail in the final segmentation. This result can be compared against a traditional Otsu’s threshold method shown in plot (f). Although the results from Figure 7e are promising, we implemented a final refinement process based on the so-called GF feathering filtering, in which Algorithm 2 is used with a carefully selected radius r and regularization parameter $ϵ$. The GF feathering filtering enables a faster implementation of Algorithm 2 by avoiding the $O\left(N^2\right)$ computational restriction of the convolution.

The quality of the proposed segmentation method was tested and compared against the binary mask segmentation of three well-known methods: (i) k-means [49], (ii) mean-shift [50], and (iii) manual threshold over the HSV color representation [51]. The acronyms used in the comparison results are listed in

Table 4. Acronyms for the segmentation methods tested.

Type	Segmentation Algorithm	Acronym
manual	HSV threshold	MCS
manual	Manual GC No refinement	GC
auto	kmeans (K = 2)	KM2
auto	meanshift (BW = 2)	MS2
auto	meanshift (BW = 4)	MS4
auto	meanshift (BW = 8)	MS8
auto	meanshift (BW = 16)	MS16
auto	GrabCut +Fixed $T\_B$ + GF	GCFauto

. By applying the fully automatic segmentation method described in Algorithm 2, we used the precision, recall, and the F1-score to measure the performance of the proposed segmentation listed as GCFauto in Table 4. Figure 8 and

Table 5. Mean Results for precision, recall, and F1-score for 10 images.

Mean	MCS	GC	KM2	MS2	MS4	MS8	MS16	GCFauto
precision	0.959	0.980	0.774	0.802	0.803	0.826	0.773	0.965
Recall	0.964	0.978	0.743	0.811	0.809	0.852	0.777	0.983
F1-Score	0.967	0.981	0.768	0.791	0.791	0.838	0.791	0.978

show the results.

For image pre-processing, the automatic GrabCut method in Algorithm 2 with the guided-filter refinement enabled precise plot segmentation of multispectral imagery with richer detail of the canopy structure and proper elimination of the background cluster. In this regard, Figure 9 shows the final segmentation results. In plot (a), the segmented image is achieved by combining the multispectral image shown in plot (b) and the RGB images from plot (c). In plot (d), the values for the radius and $ϵ$ depend on the image size. For this application, those values were experimentally determined as $r=60$ and $ϵ={0.001}^2$. Without this segmentation method, significant image corrections would be needed.

3.3. Nitrogen Estimations

In the following, we present a comprehensive comparison among multi-variable linear regressions (MLR), support vector machines (SVM), and artificial neural networks (NN) for the estimation of the canopy N. All these models were trained using the ground-truth (cf. Table 2) containing the direct measurements of leaf chlorophyll based on SPAD readings. In this regard, the estimation results are all given in SPAD scale, in which the linear relationship with nitrogen is defined as: $N=0.079\left(SPAD\right)-0.154$, according to the work reported in [37].

3.3.1. MLR Models

Figure 10 shows the estimation results using the MLR. The samples-axis corresponds to the aerial imagery used for the estimation of nitrogen thought the phenological cycle. As previously mentioned, direct SPAD measurements were conducted for selected crop plots in order to assemble the ground-truth dataset. Given that, our system selects those aerial samples matching with the GPS coordinates of the ground measurements. Overall, the MLR achieved an average N correlations of $0.93$ (V), $0.89$ (V), and $0.82$ (Ri). The data was filtered using a Moving Average Filter to reduce signal fluctuations and smooth the estimation result.

Table 6. Numerical data obtained from the multi-variable linear regressions (MLR).

Crop Stage	Correlation	RMSE
Vegetative	0.929	3.914
Reproductive	0.890	4.356
Ripening	0.822	2.318

details the numerical data. Since the coefficients for the linear regressions were found and calibrated for each stage independently, we found strong linear relationships in the vegetative stage between the VIs and the N accumulation. Through the ripening stage, the yellow color becomes predominant and parcels cannot be distinguished with ease (as shown in Figure 2c). This changes in the canopy color and dense biomass accumulation tend to saturate the linear relationships between the VIs and the canopy reflectances.

3.3.2. SVM Models

We first tried to find a single SVM model for the entire crop phenological cycle; however, the canopy reflectances and the VIs changed drastically from vegetative through ripening. In this regard, an SVM model was defined for each crop stage independently, as shown in

Table 7. Numerical data obtained from support vector machine (SVM) with several kernels.

Kernel	Crop Stage	RMSE	Correlation	${\mathit{R}}^2$
	Vegetative	2.34	0.96	0.927
Quadratic	Reproductive	2.33	0.94	0.892
	Ripening	2.04	0.87	0.770
	Vegetative	3.75	0.90	0.827
Cubic	Reproductive	12.28	0.34	0.119
	Ripening	4.39	0.57	0.326
	Vegetative	4.80	0.86	0.744
Linear	Reproductive	5.88	0.81	0.656
	Ripening	2.38	0.83	0.689
	Vegetative	5.42	0.84	0.708
Medium Gaussian	Reproductive	3.03	0.92	0.887
	Ripening	2.35	0.84	0.713
	Vegetative	5.87	0.79	0.634
Coarse Gaussian	Reproductive	4.60	0.87	0.726
	Ripening	2.97	0.81	0.660
	Vegetative	6.07	0.76	0.579
Fine Gaussian	Reproductive	4.02	0.91	0.882
	Ripening	2.16	0.87	0.761

. In addition, Figure 11 details the estimation results for the N dynamics. Based on these results, the following SVM configurations are selected for our estimation system:

• Vegetative: Quadratic, correlation = 0.96, RMSE = 2.34, $R^2$ = 0.927, $ϵ=2$.
• Reproductive: Quadratic, correlation = 0.94, RMSE = 2.33, $R^2$ = 0.892, $ϵ=2$.
• Ripening: Quadratic, correlation = 0.87, RMSE = 2.04, $R^2$ = 0.77, $ϵ=0.5$.

3.3.3. NN Models

Neural networks (NN) with several hidden layers and optimization algorithms were tested. We highlight the results with two and five hidden layers. In Figure 12, the NN with two layers achieved a correlation = 0.964, RSME = 3.315, and $R^2$ = 0.94. Increasing up to five layers, the NN achieved a correlation = 0.986, RSME = 1.703, and $R^2$ = 0.97. Several training functions were tested for each crop stage, where the BFGS Quasi-Newton functions achieved accurate results for most of the vegetative stages, whereas the Levenberg–Marquardt function for both reproductive and ripening. The numerical data are reported in

Table 8. Numerical data obtained from the 5-layer NN with several training functions.

NN Training Function	Crop Stage	Correlation	RMSE	${\mathit{R}}^2$
	Vegetative	0.979	1.986	0.959
Bayesian regression	Reproductive	0.921	3.779	0.848
	Ripening	0.67	4.405	0.462
	Vegetative	0.985	1.687	0.971
BFGS Quasi-Newton	Reproductive	0.938	3.203	0.880
	Ripening	0.851	2.149	0.724
	Vegetative	0.983	1.855	0.966
Levenberg-Marquardt	Reproductive	0.9442	3.007	0.891
	Ripening	0.890	1.835	0.792
	Vegetative	0.952	2.952	0.907
Scaled Conjugate Gradient	Reproductive	0.925	3.446	0.856
	Ripening	0.862	2.115	0.742

. Based on the results, the following NN configurations are selected for our estimation system:

• Vegetative: BFGS Quasi-Newton, correlation = 0.986, RMSE = 1.643, $R^2$ = 0.972.
• Reproductive: Levenberg-Marquardt, correlation = 0.944, RMSE = 3, $R^2$ = 0.891.
• Ripening: Levenberg-Marquardt, correlation = 0.890, RMSE = 1.835, $R^2$ = 0.792.

3.3.4. Machine-Learning Comparative Results

Figure 13 shows overall N estimation results obtained for each machine learning model. In general, we demonstrated strong correlations between the canopy N and the corresponding vegetative indices. In the case of MLR models, the coefficients for the regressions were determined and independently calibrated for each crop stage. We encountered that the N dynamics have strong linear dependencies with MSAVI, GNDVI, and NDVI; concretely, during vegetative and reproductive stages.

Table 9. Overall numerical results for canopy N estimations.

	N Estimations			N SPAD Measurements			Correlation
Methods	V	R	Ri	V	R	Ri	V	R	Ri
MLR	36.923	30.766	42.669				0.935	0.890	0.82
SVM	35.7472	29.2902	42.5226	35.7906	29.6338	42.4207	0.9699	0.9467	0.8770
NN	36.173	30.495	43.2648				0.986	0.9442	0.890

reports the overall numerical data.

From the ROC curve reported in Figure 13d, the accuracy (ACC) was calculated for each ML model. Neural networks achieved an average ACC = 0.85, improving the estimations of canopy N over the other ML models. In fact, by comparing the N-to-SPAD correlations reported on Table 9, the correlation metric was improved over each crop stage. Table 9 compares the mean N estimations achieved by each ML method against the mean SPAD-based N readings measured at ground level. Mean results are presented for each crop stage: vegetative (V), reproductive (R), and ripening (Ri). The last columns in Table 9 report the mean linear correlations between estimations and measurements.

By comparing the neural networks (NN) against the multi-variable linear regressions (MLR) and support vector machines (SVM), the former estimator achieves higher correlations partly due to the reliability and size of the assembled ground-truth database for training, i.e., reliable N datasets (direct N-leaf measurements) with sufficient variations of the crops during the growth stage (see the experimental protocol in Section 2). In addition, some VIs, such as the simple ratio (SR), tend to saturate as long the the crop grows, but other vegetative indices that behave better with higher biomass and N concentration compensate the effect of the saturated ones. The nonlinear combinations of these VIs enable the NN to achieve accurate estimations.

Furthermore, we tested the effects of the UAV control architecture proposed in Section 2.3, for improving the N estimations by means of precise position tracking during flight. Figure 14 presents the results. Plot (b) shows the X-Y tracking trajectory for each flight control scenario. The UAV was flying at a constant altitude reference of 20 m over the crop at a maximum speed of 1.5 m/s. The black boxes represent the crop plots with ground-level markers where direct N measurements were taken via SPAD. The green dotted circles show the points where the UAV stops and hovers during 4 s to capture canopy imagery.

The paths followed by the UAV are shown for both controllers; one driven by the PID that comes standard with the UAV autopilot, and the other driven by the proposed BS+DAF. Likewise, the star markers represent the GPS points of the images taken while hovering (black color for the PID and red for the BS + DAF). As observed, the attitude stabilization of the UAV seems to be more affected with the PID, since the aerial samples within each hovering area (black star markers) are spread over different positions. Instead, the aerial samples obtained with the BS + DAF (red star markers) tend to be more concentrated within the crop plot of interest. The goal is actually capturing most images in the GPS coordinates that accurately matches with the GPS position of the white markers at ground-level, as depicted in Figure 14a.

Finally, Figure 14c presents the N estimation results achieved under both control flight scenarios. Larger fluctuations, noise, and estimation errors are presented when the ML methods are trained and tested with the imagery dataset captured by the PID-driven UAV. This is mainly due to the imprecise X-Y positioning of the aerial vehicle during hovering (caused by improper attitude stabilization), in which several images capture crop areas that do not necessarily match with the ground-level markers, as detailed in Figure 14b. On the other hand, the proposed BS+DAF controller demonstrated being suitable for the proper regulation of the angular motions of the UAV by compensating external disturbances faster and achieving accurate X-Y positioning. As a result, The N estimations correlate higher with the ground-truth. As shown in the histograms in Figure 14d, the N-to-SPAD correlations grouped higher for the experiments driven by the proposed BS+DAF controller.

4. Discussion

The results presented in this work demonstrate a significant improvement in leaf nitrogen prediction in rice crops from images obtained via UAVs, validated with high SPAD correlations. These results are thanks to a combination of technological contributions, namely: (i) a novel UAV stabilization control scheme named BS+DAF that enables precise multispectral aerial image acquisition and registration with ground-truth fiducials, even in the presence of strong wind perturbations, (ii) an automated GrabCut image segmentation method, which leads to finer detail of the plant’s canopy in the RGN space, as detailed in Figure 5b and Figure 9d, and (iii) the successful application of machine learning methods trained with the vegetation indices (VI) extracted from the segmented multispectral imagery. Figure 5c confirms the increased reliability and accuracy per VI by calculating the data variance at each crop stage.

In a similar manner, Table 9 summarizes the higher nitrogen-to-SPAD correlations achieved with the implemented neural networks, and confirms the significant nonlinear relationships between leaf spectral reflectances (VI) and chlorophyll-based estimated N concentrations. These N estimation results are comparable to other state-of-the-art results presented in the scientific literature. In [52], Figure 7b shows N estimations in rice using linear regression models, for which the authors report: $RMSE=0.212$, correlation $r=0.89$, and $R^2=0.803$. Our results (tabulated in Table 8 and Table 9) reach mean correlations for the vegetative stage ($r=0.986$), reproductive ($r=0.9442$), and ripening stage ($r=0.89$). In [53], Table 6 and Figure 5 present N-status estimations in rice using vegetation indices calculated with aerial UAV imagery. Authors compared several regression models, with a highest $R^2$ of $0.74$ using the LDM method. In our case, Table 8 lists higher $R^2$ values from the Levenberg–Marquardt method.

Our approach contributes novel solutions to commercial UAV-based phenotyping technologies, particularly to image-based remote sensing applications that adopt photogrammetric geometric image post-processing methods for image correction, and enables crop/plot analysis in real-time.

One important drawback of our solution was the estimation of canopy N for crops with higher biomass. Counter to the expectation of finding a higher N correlation as a function of crop growth, the correlation results decreased in the ripening stage, as depicted in Figure 13e and from the numerical data presented in Table 9. We attribute this to the disproportionately high number of senescent leaves (yellow color with a bandwidth of 570–590 nm) that do not properly match with the VI dominant wavelengths (see Table 3). Consequently, we suggest introducing plant senescence reflectance as an added variable for the discovery of novel vegetative indices to enhance our system.

5. Conclusions

This paper presents an integrated UAV system with non-invasive image-based methods for the estimation and monitoring of the N dynamics in rice crops. Several challenges were addressed, associated with image segmentation and UAV navigation control, to achieve reliable machine learning training and highly correlated results to SPAD. The proposed BS+DAF attitude controller led to precise UAV way-point tracking, with position errors below $2\%$. This was key to capture aerial multispectral imagery during hovering that accurately matched with the GPS positions of the SPAD measurements at ground-level. This reduces the need for computationally expensive photogrammetry methods.

Higher N correlations were achieved with neural networks: $0.98$ in the vegetative stage, $0.94$ in reproductive, and $0.89$ in ripening, with an ROC accuracy of $0.85$. This is a promising result towards the autonomous estimation of rice canopy nitrogen in real-time. With the advent of small AI-dedicated systems on chip (SoC) and the powerful computing capabilities of our UAV system, we expect upcoming work to achieve real-time computation of the machine learning methods presented in this work.

Several challenges still remain for improving the N estimations, especially when the crop biomass is high. Rice crops commonly have several mixed varieties per small plot areas. Therefore, besides the assessment of the N dynamics as a function of the crop stages, it is also crucial to identify and classify the plants according to their genotype, given that nitrogen is highly sensitive to plant oxygen consumption, root growth, among other variables. To this purpose, the aerial identification of different plant varieties requires sensing spatial resolutions below 30 cm. In addition, high-dimensional data clustering is needed for the extraction of relevant features according to the plant variety.