Accuracy Comparison of Estimation on Cotton Leaf and Plant Nitrogen Content Based on UAV Digital Image under Different Nutrition Treatments

Abstract

The rapid, accurate estimation of leaf nitrogen content (LNC) and plant nitrogen content (PNC) in cotton in a non-destructive way is of great significance to the nutrient management of cotton fields. The RGB images of cotton fields in Shihezi (China) were obtained by using a low-cost unmanned aerial vehicle (UAV) with a visible-light digital camera. Combined with the data of LNC and PNC in different growth stages, the correlation between N content and visible light vegetation indices (VIs) was analyzed, and then the Random Forest (RF), Support Vector Machine (SVM), Back Propagation Neural Network (BP), and stepwise multiple linear regression (SMLR) were used to develop N content estimation models at different growth stages. The accuracy of the estimation model was assessed by coefficient of determination (R²), root mean squared error (RMSE), and relative root mean square error (rRMSE), so as to determine the optimal estimated growth stage and the best model. The results showed that the correlation between VIs and LNC was stronger than that between PNC, and the estimation accuracy of different models decreased continuously with the development of growth stages, with higher estimation accuracy in the peak squaring stage. Among the four algorithms, the best accuracy (R² = 0.9001, RMSE = 1.2309, rRMSE = 2.46% for model establishment, and R² = 0.8782, RMSE = 1.3877, rRMSE = 2.82% for model validation) was obtained when applying RF at the peak squaring stage. The LNC model for whole growth stages could be used in the later growth stage due to its higher accuracy. The results of this study showed that there is a potential for using an affordable and non-destructive UAV-based digital system to produce predicted LNC content maps that are representative of the current field nitrogen status.

1. Introduction

Cotton is an important fiber and oil crop in the world, and its planting area ranks first among cash crops in China. To ensure high and stable cotton production, chemical fertilizer application is a necessary and key step used by local farmers, while excessive application of fertilizers leads to higher costs, leaching, and soil environmental pollution, especially nitrogen (N) fertilizers, which are widely and waywardly used in cotton growth [1]. N content, as an important agronomic index and critical factor affecting cotton growth, is normally presented for cotton growth condition and yield prediction, as well as used for cotton nutrition diagnosis [2,3]. Therefore, developing a rapid and accurate estimation approach for cotton N content is of great significance for smart nutrition management in precision agriculture and also contributes to environmental protection and farmers’ economic benefits.

The traditional method to evaluate crop N status is chemical analysis, which is a destructive and time-consuming approach. Additionally, using handheld sensors (SPAD(Japan), Greenseeker(America), Li-cor(America) et al.) can indirectly evaluate crop N status at the small scale level, which cannot comprehensively represent crop canopy in the field or at the large scale level [4,5,6]. Remote sensing technology based on an unmanned aerial vehicle (UAV) has been widely used for image collection to evaluate crop growth conditions [7,8,9,10,11,12]. Most studies used RGB, near-infrared, multispectral, hyperspectral, and thermal infrared cameras or sensors based on UAVs for crop monitoring at a large scale [8,9,10,11,12,13,14]. Due to its affordability, flexibility, and high efficiency, as well as its band alignment, the RGB camera is mostly used in crop monitoring compared to multispectral and hyperspectral sensors, which indicates that the RGB camera based on UAVs has great potential for crop monitoring [15,16].

Leaf N content (LNC) or plant N content (PNC) is related to the capability of photosynthesis, which is highly affected by N fertilizer applications and crop yield; thus, many studies have been conducted to monitor the N nutrition status in leaves or plants of different crops by using UAV [2,13,17,18,19]. Jiang et al. [17] have shown that UAV-based RGB cameras have the capability to derive optimum vegetation indices (VIs) for LNC monitoring in winter wheat, and the study constructed a new VI (true color VI) from RGB images that can better mitigate saturation than other VIs. Lebourgeois et al. [16] used RGB and NIR-G-B sensors mounted on UAVs to detect N status and found the best correlation of LNC in sugarcane with a broadband version of the simple ratio pigment index (SRPIb) (R² = 0.7) among all indices examined. Further, Schirrmann et al. [20] found the red and green ratio from UAV RGB imagery correlated well (R² = 0.68) with PNC at the heading stage of winter wheat.

The aforementioned studies proved that UAV RGB imagery is a good tool for crop N nutrition status detection, but whether the cotton N content of leaves or plants at different growth stages could be detected by RGB sensors based on UAV still remains unclear. Most studies on the estimation of cotton N status used hyperspectral or multispectral sensors, which contain a large amount of band information and data for analysis, and extracting useful data from the complex data to construct a prediction model is always a tough technical challenge [21,22,23,24]. Moreover, those studies on the estimation of cotton N status using an RGB camera sensor were mostly carried out under different N treatments without considering the other fertilizer elements [25,26]. As the soil environment is not composed of simple major elements, it contains N, phosphorus, and potassium, three major elements, and different fertigation managements affect cotton N absorption and utilization, which could affect the plant N content. However, the cotton LNC or PNC estimation by using digital cameras under different nutrition conditions is poorly understood, and few studies have used digital camera-extracted data combined with machine learning technology to estimate the LNC or PNC of cotton at different growth stages.

Therefore, the objective of this work is to estimate the N content of drip-irrigated cotton in leaves and at plant level under combined N, phosphorus, and potassium treatments. To achieve this objective, thirty VIs extracted from visible digital images of low-cost UAVs, combined with three machine learning techniques and stepwise multiple linear regression (SMLR), were used to evaluate their performance in estimating N content in leaves and plants. The results of this study will provide a reference for precision agro-ecological applications and management of drip irrigation cotton fields in Xinjiang.

2. Materials and Methods

2.1. Experimental Design

Field experiments were carried out at the experiment farm of Shihezi University, Shihezi, Xinjiang, China (44°18′ N, 86°02′ E) in 2018, 2019, and 2020. A local cotton cultivar, Lumianyan 24, currently used in local production, was grown in the cotton field. The average soil fertility during 2018–2020 in the experimental site soil was alkaline and contained N 145.47 mg kg⁻¹, Olsen P 36.18 mg kg⁻¹, K₂O 67.30 mg kg⁻¹, and organic matter 31.50 g kg⁻¹ at a pH of 8.15.

These experiments were designed using a randomized complete block with three replicates. Four N application rates (506, 402.5, 299, and 195.5 kg ha⁻¹ designated as N1, N2, N3, and N4, respectively) and four different PK-Ms were designated as PK-M1 (the proportion of P and K application during squaring and bloom-bolling stages were 100%:0%), PK-M2 (the proportion of P and K application during squaring and bloom-bolling stages were 25%:75%), PK-M3 (the proportion of P and K application during squaring and bloom-bolling stages were 50%:50%), and PK-M4 (the proportion of P and K application during squaring and bloom-bolling stages were 75%:25%) (

Table 1. The experimental design in this study.

N Treatments	PK-Ms	Proportion for Squaring Stage	Proportion for Bloom-Bolling Stages
N1 (506 kg ha−1)	PK-M1	100%	0%
	PK-M2	25%	75%
	PK-M3	50%	50%
	PK-M4	75%	25%
N2 (402.5 kg ha−1)	PK-M1	100%	0%
	PK-M2	25%	75%
	PK-M3	50%	50%
	PK-M4	75%	25%
N3 (299 kg ha−1)	PK-M1	100%	0%
	PK-M2	25%	75%
	PK-M3	50%	50%
	PK-M4	75%	25%
N4 (195.5 kg ha−1)	PK-M1	100%	0%
	PK-M2	25%	75%
	PK-M3	50%	50%
	PK-M4	75%	25%

). Phosphorus and potassium fertilizers were applied to each plot, with P₂O₅ and K₂O at 108 kg ha⁻¹ and 97.2 kg ha⁻¹, respectively. There were 48 plots in the experiments, and the size of each plot was 15 m × 2.25 m.

A plot seeder was used to sow seeds with a density of 11 plants per square meter in three seasons. The sowing dates in the three seasons were 21 April 2018, 24 April 2019, and 23 April 2020, respectively. The topping was conducted on 15 July 2018, 7 July 2019, and 10 July 2020. All fertilizer is applied with water every 7–10 days during the growth period, for a total of eight times.

2.2. Data Collection

In this study, 11 ground control points were marked on the ground in advance, and the coordinates of each control point were measured with GPS. The RGB camera from the DJI Mavic Pro (

Table 2. Main parameters of UAV editions for the cotton seasons.

Edition of UAV	Camera Model	Image Size	Maxium Flight Time/min	Image Format	Weight/g
DJI Mavic pro	FC220	4000 × 3000	27	JPEG	743

and Figure 1) was used to collect data at each key growth stage of cotton. We chose the environmental conditions of stable light intensity, clear, cloudless weather, and wind speed less than level 4 to carry out continuous monitoring in the study area. Each flight activity was conducted from 12:00 to 14:00 local time. DJI GS PRO software was used to control the UAV, and each time according to the pre-defined flight plan with approximately 80% forward overlapping and 60% side overlapping, it acquired about 350 images (a spatial resolution of 1.3 cm). The frequency of image acquisition was 1 frame per 5 s, and the images were saved in JPEG format. Furthermore, the setting of the camera’s aperture was f/5. The UAV was set to automatic flying mode, flying at an altitude of 10 m and a speed of 2 m per second. Each flight activity took about 12 min. The same flight routes and camera settings were used in each growing season. According to data processing from Lu et al. (2019) [27], the acquired UAV digital images were used for image alignment, geographic reference, mosaicking, generation of dense point clouds, orthorectified image generation, and ground control point correction by Agisoft PhotoScan Professional software (

Table 3. Processing steps with corresponding parameter settings in Agisoft Photoscan software for generation of orthophotos from UAV imagery.

Task	Parameter Setup
Aligning image	Accuracy: high; Pair selection: generic; Key points: 40,000; Tie points: 4000
Building mesh	Surface type: height field; Source data: dense cloud; Face count: high
Positioning a guided marker	Manual positioning of markers on the even 11 GCPs for all the photos
Optimizing cameras	Default settings
Building dense point clouds	Quality: high; Depth filtering: mild
Building DEM	Surface: mesh; Other parameters: default
Building orthomosaic	Surface: mesh; Other parameters: default

) and ENVI. Firstly, the overlapping images were aligned automatically by using a feature point matching algorithm in software, and then all ground control points were used to georeference each image, and the internal parameters were estimated based on the image alignment and the ground control points. Secondly, the estimated parameters were used to compensate for a linear model misalignment while georeferencing the model, and we chose median filtering to build a dense point cloud. Lastly, the orthophotos were generated. Before the R, G, and B values were extracted from each plot in each growth stage, the image segmentation was processed by the Super-Green method (2 g-r-b) to obtain the Super-Green images, and then the image was processed in grayscale and enhanced by median filtering (Figure 2). The gray-level image was transformed into a binary image by fixed threshold segmentation (the value was 0.40). The color image with a black background was obtained by color filling, and the R, G, and B values in the three bands were separately set and extracted from the plot by using MATLAB R2018a (Mathworks Inc., Natick, MA, USA). After we extracted the R, G, and B values, we used the average data for further data analysis.

Plants from 48 plots were randomly sampled within a day before the UAV campaigns in each growth stage (

Table 4. Summary of field campaigns for the cotton seasons.

Year	Date of UAV Flights	Date of Field Sampling	Growth Stage
2018	24 June	23 June	Peak squaring
	13 July	12 July	Peak flowering
	2 August	2 August	Early peak boll-forming
	12 August	12 August	Late peak boll-forming
	23 August	22 August	Boll-opening
2019	25 June	24 June	Peak squaring
	19 July	19 July	Peak flowering
	6 August	5 August	Early peak boll-forming
	20 August	19 August	Late peak boll-forming
	23 September	23 September	Boll-opening
2020	14 June	14 June	Peak squaring
	4 July	3 July	Peak flowering
	27 July	26 July	Early peak boll-forming
	2 August	2 August	Late peak boll-forming
	17 September	16 September	Boll-opening

). Three representative cotton plants in each plot were selected as one sample, put into paper bags, and brought back to the laboratory for follow-up treatment. The cotton plant was divided into leaves, stems, and reproductive organs and oven-dried at 105 °C for 30 min and afterwards at 80 °C until the dry matter weight of each organ was reached.

The oven-dried samples of different organs used to determine N content (Kjeldahl method) [28]. Samples of 0.1 g of dried and ground samples were digested using a mixture of H₂O₂ and H₂SO₄ until the solution became clear and transparent, and the N content was determined using a continuous-flow auto-analyzer (BRAN + LUEBBE AA3; Germany). The formula for calculating plant N content is as follows:

Plant N content (PNC, %) = (Leaf N content (%) × Leaf dry matter weight + Stem dry N content (%) × Stem dry matter weight + Reproductive organs N content (%) × Reproductive organs dry matter weight)/Aboveground dry matter weight.

2.3. Vegetation Indices

Based on the preprocessed UAV digital image, the 48 RGB values in each growth stage are extracted by MATLAB. According to the previous research results on the relationship between N and visible light band VIs, 30 VIs were selected to estimate the LNC and PNC content of cotton. The visible light VIs selected in this study are listed in

Table 5. Summary of vegetation indices derived from aerial orthophotos for the estimation of leaf nitrogen content and plant nitrogen content in cotton.

Variable Coding	Index	Formula	References
V1	R	DN values of the R band	-
V2	G	DN values of the G band	-
V3	B	DN values of the B band	-
V4	r	R/(R + G + B)	[29]
V5	g	G/(R + G + B)	[29]
V6	b	B/(R + G + B)	[29]
V7	GRRI	G/R	[30]
V8	GBRI	G/B	[30]
V9	RBRI	R/B	[30]
V10	ARG	(R + G + B)/3	-
V11	GRVI	(G − R)/(G + R)	[31]
V12	NDI	(G − R)/(G + R)	[32]
V13	WI	(G − B)/ABS(R − G)	[29]
V14	IKAW	(R − B)/(R + B)	[33]
V15	GLI	(2 × G − R − B)/(2 × G + R + B)	[34]
V16	VARI	(G − R)/(G + R − B)	[35]
V17	ExR	1.4 × r − g	[36]
V18	ExG	2 × g − r − b	[36]
V19	ExB	1.4 × b − g	[32]
V20	ExGR	ExG − ExR	[37]
V21	VEG	g/(ra × b(1 − a)) where a = 0.667	[38]
V22	CIVE	0.441 × R − 0.881 × G + 0.385 × B + 18.78745	[39]
V23	MGRVI	G2 − R2/G2 + R2	[40]
V24	RGBVI	(G2 − B × R)/(G2 + B × R)	[40]
V25	RGRI	R/G	[41]
V26	r/b	r/b	-
V27	r − b	r − b	-
V28	r + b	r + b	-
V29	g − b	g − b	-
V30	(r − g − b)/(r + g)	(r − g − b)/(r + g)	-

.

2.4. Regression Techniques

In this study, three machine learning algorithms—Random Forest (RF), Support Vector Machine (SVM), Back Propagation Neural Network (BP) and stepwise multiple linear regression (SMLR)—were conducted by Data Processing System software.

The two important parameters (mtry and ntree) in the RF were set to 10 and 1000; the Support Vector Machine chose the radial basis function (RBF) to solve the nonlinear relationship in the data, and the cost value was set to 1. The BP neural network adopted a three-layer structure: the input layer, the hidden layer, and the output layer are respectively set up with 30-21-1, the number of iterations is 1000, and the training rate is 0.1.

2.5. Accuracy Assessment

From the samples collected in each growing season, we pooled the data from 3 years and all growing conditions and stages to form a dataset. We used the hold-out method to split the dataset into two parts, with 2/3 for model calibration and 1/3 for model validation.

In this study, the performance of each regression model was compared by three indexes: determination coefficient (R²), root mean square error (RMSE), and relative root mean square error (rRMSE). R² is used to represent the fitting effect between the predicted value and the measured value, and the closer the value is to 1, the higher the fitting accuracy of the model. The root mean square error reflects the discreteness between the predicted value and the measured value, so the smaller the value, the better the fitting accuracy of the estimation model is. The larger the R² of the inversion model and its verification, and the smaller the RMSE and rRMSE, the better the estimation ability of the model. The calculation formula is as follows:

$${\mathrm{R}}^2=\frac{{\displaystyle {\sum }_{\mathrm{i}=1}^{\mathrm{n}}({\mathrm{X}}_{\mathrm{i}}-{\overline{\mathrm{X})}}^2{({\mathrm{Y}}_{\mathrm{i}}-\overline{\mathrm{Y}})}^2}}{\mathrm{n}{\displaystyle {\sum }_{\mathrm{i}=1}^{\mathrm{n}}{({\mathrm{X}}_{\mathrm{i}}-\overline{\mathrm{X}})}^2{\displaystyle {\sum }_{\mathrm{i}=1}^{\mathrm{n}}{({\mathrm{Y}}_{\mathrm{i}}-\overline{\mathrm{Y}})}^2}}}$$

(1)

$$\mathrm{RMSE}=\sqrt{\frac{{\displaystyle {\sum }_{\mathrm{i}=1}^{\mathrm{n}}{({\mathrm{Y}}_{\mathrm{i}}-{\mathrm{X}}_{\mathrm{i}})}^2}}{\mathrm{n}}}$$

(2)

$$\mathrm{rRMSE}(100\%)=\frac{\mathrm{RMSE}}{\overline{\mathrm{X}}}\times 100$$

(3)

where ${\mathrm{X}}_{\mathrm{i}}$, $\overline{\mathrm{X}}$, ${\mathrm{Y}}_{\mathrm{i}}$, and $\overline{\mathrm{Y}}$ represent the measured value, the mean value, the average value of the predicted value, and the average value of the predicted value, respectively, and n represents the sample number of the developed model.

3. Result

3.1. Variability in the Measured N Content and RGB

Figure 3 shows the overall trend of LNC and PNC status from the five growth stage data sets during 2018–2020. The trend of LNC and PNC decreased along with growth stages, and peak points were at the peaking squaring stage. The LNC varied from 16.80 g kg⁻¹ to 59.11 g kg⁻¹ across growth stages. Average LNC decreased from 49.75 g kg⁻¹ at peak squaring to 27.20 g kg⁻¹ at boll-opening stages. The PNC varied from 6.08 g kg⁻¹ to 49.89 g kg⁻¹ across growth stages. Average PNC decreased from 40.21 g kg⁻¹ at peak squaring to 17.48 g kg⁻¹ at boll-opening stages. The coefficient of variation in LNC and PNC increased from peak squaring stages to boll-opening stages (

Table 6. Basic statistics of the nitrogen content in cotton at different growth stages.

Parameters	Stage	Min (g/kg)	Max (g/kg)	Mean (g/kg)	CV %
LNC	Peak squaring	40.73	59.11	49.75 ± 3.84	7.73
	Peak flowering	31.11	52.46	44.54 ± 4.07	9.13
	Early peak boll-forming	28.41	48.35	38.65 ± 4.31	11.15
	Late peak boll-forming	21.74	41.46	33.41 ± 3.84	11.49
	Boll-opening	16.80	39.58	27.20 ± 4.28	15.72
PNC	Peak squaring	32.06	49.89	40.21 ± 3.51	8.73
	Peak flowering	26.57	40.37	34.57 ± 3.31	9.57
	Early peak boll-forming	19.82	35.50	28.32 ± 3.36	11.87
	Late peak boll-forming	14.08	32.37	23.60 ± 3.85	16.30
	Boll-opening	6.08	26.83	17.48 ± 4.00	22.91

). Based on the analysis of variance, the year, N treatment, and PK-Ms significantly affected LNC and PNC, and there was an interaction between N and PK-Ms on LNC and PNC during 2018–2020 (

Table 7. Analysis of variance for leaf nitrogen content and plant content at nitrogen treatments and phosphorus-potassium managements (PK-Ms) during 2018–2020.

Treatment	LNC	PNC
Year	**	**
N	**	**
PK-Ms	**	**
Year × N	**	**
Year × PK-Ms	**	**
N × PK-Ms	**	*
Year × N × PK-Ms	ns	ns

Note: ** means significant at the level of 0.01’ ; * means significant at the level of 0.01; ns means not significant.

).

Figure 4 and

Table 8. Basic statistics of the R, G, and B values in cotton at different growth stages.

Band	Stage	Min	Max	Mean	CV (%)
R	Peak squaring	77.20	114.99	94.74 ± 7.05	7.45
	Peak flowering	73.53	126.27	102.29 ± 9.27	9.06
	Early peak boll-forming	97.24	134.89	113.94 ± 8.08	7.09
	Late peak boll-forming	105.31	143.29	124.56 ± 9.06	7.27
	Boll-opening	113.44	140.10	125.71 ± 5.64	4.48
G	Peak squaring	89.85	166.29	121.49 ± 14.29	11.76
	Peak flowering	72.79	170.78	120.81 ± 17.68	14.63
	Early peak boll-forming	98.13	179.96	132.6 ± 15.22	11.48
	Late peak boll-forming	110.54	183.46	141.88 ± 13.92	9.81
	Boll-opening	114.71	144.92	126.63 ± 6.01	4.75
B	Peak squaring	70.33	104.89	86.50 ± 6.39	7.39
	Peak flowering	65.70	114.86	92.08 ± 8.30	9.02
	Early peak boll-forming	87.82	124.12	103.66 ± 7.35	7.09
	Late peak boll-forming	97.76	132.78	115.08 ± 8.83	7.68
	Boll-opening	112.83	140.96	124.97 ± 5.63	4.50

show the overall trend of R, G, and B status from the five growth stage data sets during 2018–2020. The trend of R value increased first and then decreased along with growth stages, except in 2018. The R value varied from 73.53 to 143.29 across growth stages. The average R value increased from 94.74 at peak squaring to 125.71 at boll-opening stages. The trend of G value increased first and then decreased at the boll-opening stage, and the G value varied from 72.79 to 183.46 across growth stages. The lowest average G value was 120.81 at the peak flowering stage, and the highest G value was 141.88 at the late peak boll-forming stage. The trend of the B value increased along with growth stages, and the B value varied from 65.70 to 140.96 across growth stages. The lowest average B value was 86.50 at the peak squaring stage, and the highest average B value was 124.97 at the boll-opening stages.

3.2. Correlation between UAV-Derived Variables

Based on the UAV RGB band information, thirty VIs were selected, referring to the previous research results and the relationship between N and visible light band VIs (Figure 5). For the LNC, twenty-three of the thirty Vis showed positive or negative correlations. GBRI, bn, ExB, g-b, INT, and r-b were the most strongly and weakly correlated to LNC, respectively (GBRI and g-b: r = 0.80, p-value < 0.01; bn and ExB: r = −0.80, p-value < 0.01; INT: r = 0.05, p-value > 0.05; r-b: r = −0.05, p-value > 0.05). For the PNC, twenty eight of the thirty VIs showed positive or negative correlations. The highest correlations with PNC were found for bn (r = 0.79, p-value < 0.01), with the lowest for r-b (r = 0.01, p-value > 0.05). Generally, the correlations of VIs with PNC were weaker than those with LNC. Among the twenty-three of thirty VIs in LNC and twenty-eight of thirty VIs in PNC, the input parameters for model development were that the correlation coefficient between VIs and LNC (PNC) was greater than 0.50.

3.3. Comparison of N Content Estimation Performance with Machine Learning Techniques and SMLR

Table 9. Cotton nitrogen content estimates derived using UAV imagery and vegetation indices combined with three machine learning techniques and stepwise multiple linear regression.

Nitrogen Content	Stage	RF			SVM			BP			SMLR
		R2	RMSE	rRMSE (%)	R2	RMSE	rRMSE (%)	R2	RMSE	rRMSE (%)	R2	RMSE	rRMSE (%)
LNC	Peak squaring	0.9001	1.2309	2.46	0.8703	1.3883	2.78	0.7566	1.9088	3.82	0.8664	1.3973	2.79
	Peak flowering	0.7967	1.8650	4.16	0.7261	2.1954	4.90	0.7123	2.2181	4.95	0.7638	1.9468	4.34
	Early peak boll-forming	0.7430	2.2689	5.84	0.7033	2.3957	6.16	0.5939	3.0686	7.89	0.7119	2.3228	5.98
	Late peak boll-forming	0.6785	2.3477	6.98	0.6717	2.4612	7.32	0.6392	2.3804	7.08	0.6744	2.1797	6.48
	Boll-opening	0.4964	3.2866	11.95	0.5191	3.2787	11.92	0.4946	3.5326	12.85	0.4927	9.1517	33.28
	All stage	0.7249	4.8741	12.51	0.6601	5.4362	13.95	0.6401	6.2970	16.16	0.6797	5.1064	13.10
PNC	Peak squaring	0.8757	1.2324	3.06	0.8556	1.3354	3.31	0.7468	1.7663	4.38	0.8356	1.4080	3.49
	Peak flowering	0.7094	1.7874	5.16	0.6930	1.8920	5.46	0.6896	1.8553	5.36	0.6871	1.8217	5.26
	Early peak boll-forming	0.6781	2.0032	7.08	0.5758	2.1281	7.52	0.5596	2.4951	8.79	0.6147	1.9852	7.02
	Late peak boll-forming	0.6461	2.4598	10.40	0.6052	2.3304	9.85	0.4816	2.6834	11.35	0.6220	2.2843	9.66
	Boll-opening	0.4138	3.4138	19.36	0.4930	2.8577	16.21	0.4325	2.9751	16.87	0.4782	2.8305	16.05
	All stage	0.7003	5.6298	19.47	0.6381	6.6715	23.08	0.6363	5.3468	18.49	0.6648	5.5019	19.03

Note: The total number of samples was 480 for all stages and 96 for each growth stage.

lists three machine learning techniques and SMLR for the estimation of LNC and PNC over the critical growth stages. For LNC, peak-squaring growth-stage-specific models yield the best performance, with R² ≥ 0.7566, RMSE ≤ 1.9088, and rRMSE ≤ 3.82% in three machine learning techniques and SMLR. For PNC, peak-squaring growth-stage-specific models yield the best performance, with R² ≥ 0.7468, RMSE ≤ 1.7663, and rRMSE ≤ 4.38% in three machine learning techniques and SMLR. For both LNC and PNC, the growth-stage-specific models at boll-opening stages show the worst performance, with R² ≤ 0.5191, RMSE ≥ 3.2787, and rRMSE ≥ 11.92%. The RF of the three machine learning techniques could yield higher R², lower RMSE, and rRMSE values than BP, SVM, and SMLR.

The estimation model of cotton N content should not only have high accuracy but also have repeatability and universality. Therefore, the validation set data were used to verify the model constructed above; the predicted value of N content and the measured value in the field were linearly fitted, and the reliability of the selected model was shown in a 1:1 scatter plot (Figure 6). Consistently, according to R², the model performances are ordered as follows: Peak squaring > Peak flowering > Early peak boll-forming > Late peak boll-forming > Boll-opening. Peak-squaring growth-stage-specific models yield the best performance for LNC and PNC in validation. RF yielded the highest accuracy for LNC (R² = 0.8782, RMSE = 1.3877, rRMSE = 2.82%) and PNC (R² = 0.8622, RMSE = 1.3188, rRMSE = 3.30%). The BP has the lowest estimation accuracy among the three models. The growth-stage-specific models at boll-opening stages yielded the worst performance for LNC and PNC, with R² below 0.5362, RMSE above 3.2019, and rRMSE above 12.04%.

4. Discussion

The goal of this study was to predict cotton N content at leaf and plant level by using UAV imagery under different fertigation management and compare the performance of the predicted model on LNC and PNC at different growth stages. LNC and PNC, being indicators for N nutrition status, were mostly used in remote sensing technology for crop N monitoring and prediction [2,7,10,15,17,18]. As crop N status affects the chlorophyll content in plants, resulting in changes in visible light absorption and reflection, crop N nutrition diagnosis can be realized by using visible light digital image information [21,22,26,42,43]. Previous studies have shown that UAV-mounted RGB sensors have the best potential for crop N content prediction, and N content is sensitive to visible bands in the whole cotton growth stage [16,17,21]. In this study, based on the results of correlation analysis, the selected VIs have a strong correlation with cotton LNC and PNC, most of them at a significant level, which indicates that there is great potential in using digital images to estimate cotton N nutrition [16,17,20,21,22].

In this study, the UAV digital image was used to construct the estimation model of N content in leaves and plants of cotton at each key growth stage. It was found that with the development of the growth stage, the accuracy of the model constructed by digital images decreased continuously. The validation results also showed the same trend, which is similar to the results of Cao et al. [44] and Bending et al. [40], that the accuracy of VIs in estimating N status and biomass in the early growth stage is higher than that in the later growth stage. The reason for the higher accuracy in N content estimation found in crop early growth stages might be due to lower chlorophyll content. As high chlorophyll content in crops at later growth stages would easily make the saturation in visible light absorption. Further, the estimation accuracy of the model based on visible light VIs in LNC is better than that of PNC in different growth stages, which resulted from the vertical observation mode adopted in this study, and the contribution rate of canopy leaves to the spectral information of the visible band is absolutely dominant; another reason for that might be that the assembly and distribution of N in cotton is accompanied by the change of plant structure, and the variability of the PNC data set is always higher than that of LNC in the reproductive growth stage (Table 5). In addition, the accuracy of the model developed by VIs in the boll-opening stage is much lower than that in the other four stages, which is due to the senescence and shedding of leaves in the boll-opening stage, the canopy image being affected by soil spectral reflection, and the mixed information contained in the image pixels being more complex. However, in cotton production, N fertilizer is usually no longer applied in the boll-opening stage, so the low accuracy model of cotton LNC and PNC in this period has no production significance, which is similar to the results in Wang et al. [21], Ma et al. [24], and Kou et al. [26] that lower accuracy in N content was found in the boll-opening stage and better prediction in N content was found in earlier growth stages. The accuracy of the model based on the all-stage model is generally better than that in the late peak boll-forming stage, so the all-stage model can be used instead of the late peak boll-forming stage model to obtain higher accuracy, which fills the gap for N content estimation in the later growth stage of cotton compared to previous research [21,22,23,24,25].

In this study, RGB and its combined VIs are used as input variables to estimate LNC and PNC for evaluating the performance of three machine learning techniques. The results showed that the performance of RF was better than that of SVM and BP neural networks, which is similar to previous research [45,46,47,48]. RF is a popular nonlinear integrated statistical method that has a good tolerance to noise, avoids overfitting, and can solve nonlinear problems through kernel functions [46,47], so it has higher fitting accuracy than SMLR. SVM was widely used in the field of remote sensing, and it was suitable for datasets with small sample sizes and high feature dimensions [48]. As one of the typical models of artificial neural networks, the BP neural network is composed of an input layer, a hidden layer, and an output layer, which can solve complicated nonlinear mapping problems [27,49]. In this study, all VIs were calculated by RGB, so there was a strong correlation between them. Further, the RF algorithm is insensitive to the collinearity of variables [50], so it is more suitable to deal with multiple interrelated variables and does not need to test the normality and independence of variables, so it can obtain higher model accuracy. The reason why the accuracy of BP was lower than that of the other models may be that its convergence speed is slow or the learning rate is unstable in the learning process.

Some studies focused on developing a model for estimation of crop N content under different N treatments or under combined N treatments and plant densities without considering the other two major elements in the field during crop growth stages [21,22,23,24,25,26,51]. While the crop N content was affected by field fertigation management under drip-irrigated conditions, there were interaction effects among N, phosphorus, and potassium in the soil, which affected crop canopy N content, canopy reflectance, and model accuracy. In this study, the LNC and PNC were predicted based on UAV images under combined N, phosphorus, and potassium treatments, and the result was satisfactory. However, this study obtained a fixed-resolution digital image at a flight altitude of 10 m, which was relatively lower than other studies [19,25,26,27]. The flight altitude is one of the most important parameters affecting the observation accuracy, and the higher the flight altitude, the lower the spatial resolution obtained. Zhu et al. (2023) and He et al. (2022) [52,53] found that the higher flight altitude led to a larger sampling area in each pixel with fewer sampling quantities, which indicated that there was a higher possibility of decreasing the frequency of feature values and affecting the accuracy of value extraction due to increased background interference. Further, flight altitudes of 10 m or less could decrease the effect of relative humidity, air mass density, and atmospheric transmittance on data acquisition but increase the flight time and cost of UAVs [52,53,54]. Our previous studies have shown that it is a useful and efficient method to set the flight altitude at 10 m with 1.3 cm ground resolution for image acquisition and cotton detection [55,56,57], while whether the different flight altitude would affect the accuracy of the cotton growth model still needs to be further discussed.

5. Conclusions

This study presented a cotton N estimation model at leaf and plant level based on UAV images under different fertigation managements and compared the estimation accuracy of the N model between leaf and plant with three machine learning algorithms and stepwise multiple linear regression at different growth stages. The following conclusions were drawn:

The cotton N content of leaf or plant level at different growth stages could be detected by RGB sensors based on UAV, and the estimation accuracy in leaf was better than that in plant with the random forest algorithm at each growth stage, and higher accuracy was found in the peak squaring stage.

The LNC content estimation model using the random forest algorithm for whole growth stages could be used in later growth stages due to its higher accuracy.

The results of this study showed that there is potential for using color indices derived from UAV-based digital systems to produce predicted LNC content maps. In future work, texture features could be taken into consideration, as could combinations of texture and color features, to improve the accuracy of the LNC estimation model for covering more growth stages.