Estimation of Paddy Rice Nitrogen Content and Accumulation Both at Leaf and Plant Levels from UAV Hyperspectral Imagery

Abstract

Remote sensing-based mapping of crop nitrogen (N) status is beneficial for precision N management over large geographic regions. Both leaf/canopy level nitrogen content and accumulation are valuable for crop nutrient diagnosis. However, previous studies mainly focused on leaf nitrogen content (LNC) estimation. The effects of growth stages on the modeling accuracy have not been widely discussed. This study aimed to estimate different paddy rice N traits—LNC, plant nitrogen content (PNC), leaf nitrogen accumulation (LNA) and plant nitrogen accumulation (PNA)—from unmanned aerial vehicle (UAV)-based hyperspectral images. Additionally, the effects of the growth stage were evaluated. Univariate regression models on vegetation indices (VIs), the traditional multivariate calibration method, partial least squares regression (PLSR) and modern machine learning (ML) methods, including artificial neural network (ANN), random forest (RF), and support vector machine (SVM), were evaluated both over the whole growing season and in each single growth stage (including the tillering, jointing, booting and heading growth stages). The results indicate that the correlation between the four nitrogen traits and the other three biochemical traits—leaf chlorophyll content, canopy chlorophyll content and aboveground biomass—are affected by the growth stage. Within a single growth stage, the performance of selected VIs is relatively constant. For the full-growth-stage models, the performance of the VI-based models is more diverse. For the full-growth-stage models, the transformed chlorophyll absorption in the reflectance index/optimized soil-adjusted vegetation index (TCARI/OSAVI) performs best for LNC, PNC and PNA estimation, while the three band vegetation index (TBVI_Tian) performs best for LNA estimation. There are no obvious patterns regarding which method performs the best of the PLSR, ANN, RF and SVM in either the growth-stage-specific or full-growth-stage models. For the growth-stage-specific models, a lower mean relative error (MRE) and higher R² can be acquired at the tillering and jointing growth stages. The PLSR and ML methods yield obviously better estimation accuracy for the full-growth-stage models than the VI-based models. For the growth-stage-specific models, the performance of VI-based models seems optimal and cannot be obviously surpassed. These results suggest that building linear regression models on VIs for paddy rice nitrogen traits estimation is still a reasonable choice when only a single growth stage is involved. However, when multiple growth stages are involved or missing the phenology information, using PLSR or ML methods is a better option.

1. Introduction

Nitrogen (N) is one of the most important plant macronutrients and plays a significant role in the photosynthesis process. It is important for crop development, production and quality [1,2,3]. Proper N fertilization can help crops achieve maximum potential yield [4]. However, the uncertainty in the climate of the current year, previous crop management, and soil N supply cause actual N uptake to vary in both time and space. N deficiency can negatively affect photosynthesis and biomass accumulation, as well as both the quality and quantity of the crop yield [5]. Thus, considering the relatively low cost of N fertilization, overapplication of N fertilizer is a common phenomenon in agricultural practice to ensure maximum yield [6]. The overapplication of N fertilizer can result in problems such as plant stress, the overproduction of leaves, and more unproductive tillers. Additionally, it increases the N leaching risk of cropping systems.

For early growth stages, monitoring the crop N content can help in fertilization planning. For maturing growth stages, monitoring crop N content can help indicate yield quality [7]. Thus, the continuous monitoring of crop N content is of great importance. The destructive determination of N content, using chemical methods such as the Kjeldahl technique, is laborious, lengthy and costly [5]. Portable chlorophyll meters, such as SPAD- 502, can be used to quickly diagnose crop N status. These tools are still point-sampling based and not suitable for precision N management across large areas. Remote sensing, however, provides a cost-effective and nondestructive avenue for monitoring crop N status over large geographic areas [8]. Using radiometric data to map crop N conditions has been reported for many crops.

Although there are a few works trying to use radiation transformation models, such as N-PROSPECT [9] or PROPSECT-PRO [10], which are modified versions of the widely used PROSPECT model [11], for N content estimation, remote sensing-based crop N content detection is mainly conducted using empirical relationships from vegetation indices (VIs) that are sensitive to chlorophyll content based on the fact that the two variables are moderately to be highly correlated [12]. Tian et al. [13] quantified the relationship between paddy rice leaf nitrogen content (LNC) and red-edge position (REP) with different algorithms. Li et al. [14] found that the canopy chlorophyll content index (CCCI) and N planar domain index (NDPI) are stable predictors for wheat plant nitrogen content (PNC) after the heading growth stage. Patel et al. [15] evaluated the relationships of several popular broadband VIs with the PNC, plant nitrogen accumulation (PNA) and aboveground biomass (ABG) of ryegrass. Li et al. [16] evaluated several hyperspectral VIs for winter wheat PNC estimation across different years and growth stages. Fitzgerald et al. [17] suggested that the CCCI is suitable for predicting rainfall wheat PNC at the stem elongation growth stage.

There are also studies trying to build VIs that are directly related to N traits. Lee et al. [18] suggested that the first derivative of canopy spectra at 735 nm is a good indicator for the PNC of paddy rice at the panicle formation stage. Stroppiana et al. [19] concluded that the normalized difference vegetation index (NDVI) formulated with narrow band reflectance at 483 nm and 503 nm was optimal for paddy rice PNC estimation in a combined data set from panicle initiation and heading stages. They also pointed out that VIs specially designed for chlorophyll/nitrogen estimation, such as the transformed chlorophyll absorption in reflectance index (TCARI) and modified chlorophyll absorption ratio index (MCARI) did not perform well in their data set. Tian et al. [20] proposed a three band index that performed well at both the ground and space scales for paddy rice LNC estimation.

Considering that VIs typically use two to three bands, univariate regression models built on VIs are criticized because they do not fully take advantage of spectrometric data (in particular, hyperspectral data) and may be too simple to characterize the intrinsic relationships between the target variable and spectrometric data [21,22,23,24]. Partial least squares regression (PLSR) is a powerful alternative method and can perform better for biochemical trait estimation in most cases [25,26,27,28,29,30]. Specifically, Hansen and Schjoerring [24] showed that PLSR could reduce the root mean square error (RMSE) of LNC estimation by 24%, compared to those of narrowband VIs. More recently, state-of-art machine learning (ML) methods, such as artificial neural networks (ANNs), random forests (RFs) and support vector machines (SVMs), were evaluated for N trait estimation and yielded encouraging results [31,32,33,34,35,36,37,38].

Most of the previous studies focused on special narrow growth stages, such as the early season [32], stem elongation stage [17], panicle initiation stage [18,26,39], heading stage [14], and anthesis stage [40]. Additionally, the LNC is usually considered the target variable [1]. However, the N content can be expressed both per dry weight and per unit surface. Additionally, the N content can be measured at the leaf scale and at the wholeplant scale. The growth stage development is accompanied by canopy structure variation, biomass accumulation and N allocation. How the prediction strength of spectrometric data for different N traits varies across growth stages has not been fully explored. Therefore, the objective of this study is to evaluate the performance of VIs, PLSR and ML methods for paddy rice N trait estimation, with a special focus on the effects of growth stages.

2. Materials and Methods

2.1. Study Area and Experimental Setup

The experiments were conducted at the rice research field of Guangdong Academy of Agricultural Sciences (23°23′38″N, 113°25′37″E). This region is characterized by a subtropical monsoon climate. The average annual precipitation and average temperature are 1700 (mm) and 22.5 °C, respectively. The parent material of the soil in the test area is lateritic sandy loam formed by river alluvium. In this region, paddy rice is under a two-season cropping system. Our experiments were conducted during the 2020 late rice growing season. The paddy rice variety Huangruanxiuzhan was used as the test material. The paddy rice was sown in a nursery bed and manually transplanted on 9 August 2020. The transplanting density was 400 plants/plot with one plant per hill.

To acquire the necessary reference observations, a paddy rice field with nitrogen treatment was set up and used to conduct field-data collection campaigns (detailed in Figure 1). Nitrogen from 0 to 180 (kg ha⁻¹) with a 30 (kg ha⁻¹) increase was applied to plots J0 to J6 7d after transplantation (DAT). Nitrogen from 0 to 90 (kg ha⁻¹) with an  15 (kg ha⁻¹) increase was applied to plots F0 to F6 at  33 DAT. Additionally, 90 (kg ha⁻¹) nitrogen was applied to each plot from F0 to F6 at 7 DAT. The paddy rice was at the tillering and jointing growth stages at  7 and  33 DAT, respectively. Phosphates and potash fertilizers were applied as basal dressings with P₂O₅ and K₂O at 45 and 120 (kg ha⁻¹), respectively, for each plot. For both J0 to J6 and F0 to F6, each nitrogen level comprised two replicates. Thus, there were 28 plots in total. Each plot was (5 × 4.5 m) in size.

2.2. Field Data Collection

Four field-data collection campaigns were carried out at 18, 32, 47 and 62 DAT. The paddy rice was in the tillering, jointing, booting, and heading growth stages on the corresponding dates. The growth stage was determined, according to the rules described by [41].

Each plot was equally split into two sampling areas (demonstrated in Figure 1). For each field-data collection campaign, the grazeable biomass of three to five hill plants of each sampling area was randomly destructively sampled. Then, the sample plants were split into two parts: green leaves and other (senescent leaves, stems and particles). For each sample, 30 random SPAD-502 (Minolta Camera Co., Osaka, Japan) readings were recorded on green leaves. Then, all green leaves were scanned by a portable scanner (PERFECTION V39, Seiko Epson Corp., Nagano, Japan). The soil plant analysis development (SPAD) readings were converted to leaf chlorophyll content (LCC) (mg cm⁻²) with an empirical relationship (Equation (1)) proposed by [42] and averaged as the LCC of the corresponding sampling area. The leaf pixels in scanned images were extracted by a threshold method. The leaf area indices (LAIs) of each sampling area were determined by Equation (2), where ${\rho }_{greenleaf}$ is the green-leaf pixel ratio of each scanned image, $are{a}_{img}=216\times 297\times {10}^{-4}\mathrm{m}$ is the maximum scan area of the scanner, and $density=400/(5\times 4.5){\mathrm{m}}^{-2}$ is the planting density. The canopy chlorophyll content (CCC) was approximated by multiplying the LAI and LCC (Equation (3)).

$$\mathrm{LCC}=99\times SPAD/(144-SPAD)$$

(1)

$$\mathrm{LAI}=\frac{1}{density}\sum ({\rho }_{greenleaf}\times are{a}_{img})$$

(2)

$$\mathrm{CCC}=\mathrm{LAI}\times \mathrm{LCC}$$

(3)

In the second step, samples were oven-dried at 105 °C for half an hour and then at 75 °C until a constant weight was reached. Finally, the oven-dried samples of different organs were weighed and used to determine the leaf biomass, ABG, nitrogen content (Kjeldahl method) and accumulation. Both LNC and PNC were measured. Additionally, LNA and PNA were determined by multiplying LNC and PNC with leaf biomass and ABG, respectively.

Plots J0 to J6 were sampled during data collection campaigns at 18, 32 and 47 DAT, while plots F0 to F6 were sampled during data collection campaigns at 47 and 62 DAT. Thus, 28 samples were acquired for each data collection campaign, except for 47 DAT (56 samples).

Canopy reflectance data were acquired by a Cubert S185 hyperspectral camera (Cubert GmbH, Ulm, Germany) installed on a hexacopter (DJI M600Pro, SZ DJI Technology Co., Shenzhen, Guangzhou, China). The hyperspectral camera has a spectral range of 450–950 nm and a spectral resolution of 4 nm. It was installed through a three-axis gimbal stabilizer (DJI Ronin-MX, SZ DJI Technology Co., Shenzhen, Guangzhou, China). All unmanned aerial vehicle (UAV) flight missions were performed at approximately 11:00 a.m. Considering that the camera viewing angle may affect the estimation of the canopy parameters [43,44], the same view angle was ensured by the gimbal stabilizer to acquire the nadir images. The flight height was 30 m, which resulted in 16 cm spatial resolution for hyperspectral images. The hyperspectral images were mosaicked and orthographically projected for further analysis. The bands beyond 850 nm were dropped because of the low spectral quality [45]. The sample-area-wise mean spectra were calculated using regions of interest generated from sample area boundaries (demonstrated in Figure 1).

2.3. Methods

Linear regression models based on VIs, as well as PLSR and three ML models (ANN, RF and SVM) based on raw reflectance spectra, were built. Both full-growth-stage and growth-stage-specific models were built for each nitrogen trait (LNC, PNC, LNA and PNA). The general workflow is demonstrated in Figure 2. All models were built in the R environment [46]. The PLSR, ANN, RF and SVM models were built using the R packages plsmod [47], nnet [48], randomForest [49] and kernlab [50], respectively.

Twenty-three VIs were evaluated in this study (listed in

Table 1. VIs evaluated in this study.

Index	Definition	Cite
NDVI	$\frac{{\rho }_{800}-{\rho }_{670}}{{\rho }_{800}+{\rho }_{670}}$	[52]
SAVI	$(1+0.5)\frac{{\rho }_{800}-{\rho }_{670}}{{\rho }_{800}+{\rho }_{670}+0.5}$	[53]
OSAVI	$(1+0.16)\frac{{\rho }_{800}-{\rho }_{670}}{{\rho }_{800}+{\rho }_{670}+0.16}$	[54]
MSAVI	$\frac{2{\rho }_{800}+1-\sqrt{{(2{\rho }_{800}+1)}^2-8({\rho }_{800}-{\rho }_{670})}}{2}$	[55]
CIgreen	$\frac{{\rho }_{800}}{{\rho }_{560}}-1$	[56]
CIre	$\frac{{\rho }_{800}}{{\rho }_{720}}-1$	[56]
NDRE	$\frac{{\rho }_{780}-{\rho }_{720}}{{\rho }_{790}+{\rho }_{720}}$	[17]
D735	${dR/d\lambda |}_{735}$	[18]
NDVI(483,503)	$\frac{{\rho }_{503}-{\rho }_{483}}{{\rho }_{503}+{\rho }_{483}}$	[19]
TBVItian	$\frac{{\rho }_{705}}{{\rho }_{717}+{\rho }_{491}}$	[20]
mND705	$\frac{{\rho }_{750}-{\rho }_{705}}{{\rho }_{750}+{\rho }_{705}+2\times {\rho }_{445}}$	[57]
mSR705	$\frac{{\rho }_{750}-{\rho }_{445}}{{\rho }_{705}-{\rho }_{445}}$	[57]
TCARI	$3[({\rho }_{700}-{\rho }_{670})-0.2({\rho }_{700}-{\rho }_{550})({\rho }_{700}/{\rho }_{670})]$	[58]
TCARI/OSAVI	$\frac{TCARI}{OSAVI}$	[58]
TCARI2	$3[({\rho }_{750}-{\rho }_{705})-0.2({\rho }_{750}-{\rho }_{550})({\rho }_{750}/{\rho }_{705})]$	[59]
OSAVI2	$(1+0.16)\frac{{\rho }_{750}-{\rho }_{705}}{{\rho }_{750}+{\rho }_{705}+0.16}$	[59]
TCARI2/OSAVI2	$\frac{TCARI2}{OSAVI2}$	[59]
MTCI	$({\rho }_{754}-{\rho }_{709})/({\rho }_{709}-{\rho }_{681})$	[60]
REPLi	REP through linear four-point interpolation	[61]
REPLE	REP through linear extrapolation	[62]
mREIP	REP through Gaussian fit	[63]
CCCI	$\frac{NDRE-min\_NDRE}{max\_NDRE-min\_NDRE}$	[17]
NDPI	$\frac{C{I}_{re}-min\_C{I}_{re}}{max\_C{I}_{re}-min\_C{I}_{re}}$	[14]

${\rho }_{\lambda }$ stands for the reflectance at band $\lambda$ nm. If $\lambda$ does not match the S185 wavelengths, the mean reflectance of two adjacent bands is used. Considering that the $\lambda$ = 445 nm band does not fall in the S185 band ranges, the closest band, $\lambda$ = 450 nm, is used.

). These indices are commonly used for biomass, chlorophyll and nitrogen status estimation. Specifically, the CCCI [17] is a two-dimensional index based on the normalized difference red edge (NDRE) and NDVI. By capturing the planar domain idea proposed by Clarke et al. [51], the CCCI attempts to separate the soil background signal from the plant signal. The NDVI is used as a surrogate for ground cover, and the NDRE is used as a measure of nitrogen content. Conceptually, this allows the nitrogen for any ground cover to be measured. Similarly, Li et al. [14] proposed the NDPI by changing the NDRE in CCCI to CI_re.

A detailed description of the PLSR, ANN, RF and SVM methods can be found in our previous work [21]. In this work, a single-hidden-layer neural network structure was used for the ANN; a radial basis kernel function (RBF) was used for the SVM. In these methods, the latent variables (LVs) of the PLSR, the number of hidden units (units) of the ANN, the number of randomly selected variables of each tree (mtry) and the number of trees (ntree) of the RF, and the σ in the RBF and the C-constant of the regularization term in the Lagrange formulation of the SVM are hyperparameters that need to be fine tuned.

The full-growth-stage models were built on the full data set, while the growth-stage-specific models were built on the data of the corresponding growth stages. The corresponding data set was randomly split into calibration and validation data sets at a 3:1 ratio. All models were calibrated and validated on the calibration and validation data sets, respectively. The mean relative error (MRE), RMSE and coefficient of determination (R²) were calculated to evaluate the model performance (Equations (4)–(6)), where O_i and P_i are the ith observed and predicted data, respectively, and n is the size of the data pairs.

$${\mathrm{R}}^2=1-\frac{{\sum }_{i=1}^n{(O_i-P_i)}^2}{{\sum }_{i=1}^2{(O_i-P_i)}^2}$$

(4)

$$\mathrm{RMSE}=\sqrt{\frac{{\sum }_{i=1}^n{(O_i-P_i)}^2}{n}}$$

(5)

$$\mathrm{MRE}=100\%\ast \frac{1}{n}\sum _{i=1}^n\frac{|O_i-P_i|}{O_i}$$

(6)

For the PLSR, ANN, RF and SVM, the hyperparameters (LVs for the PLSR, units for the ANN, mtry for the RF and C, σ for the SVM) were determined by a grid search with a repeated (5 times) 10-fold cross-validation procedure on the calibration data set (Algorithm 1). The parameter (or parameter combination) associated with the lowest MRE_cv was considered optimal. After the key parameters were determined, the final models were calibrated on the corresponding calibration data set and evaluated on the corresponding standalone validation data set. Although no hyperparameters needed to be determined for the linear models based on VIs, the same cross-validation procedure was used to calculate the cross-validated calibration MRE_cv, RMSE_cv, R${}_{cv}^2$. The same fold split scheme across different models (of the same growth stage) was used to ensure a fair comparison. MRE_cv, RMSE_cv and R${}_{cv}^2$ were mainly used for model performance comparison.

Algorithm 1: Hyperparameter optimization and model calibration algorithm.

3. Results

3.1. Variability in the Measured Nitrogen, Chlorophyll and Aboveground Biomass

Table 2. Statistics of the observed LNC, %; PNC, %; LNA, g m⁻²; PNA, kg ha⁻¹; LCC, μg cm⁻²; CCC, μg cm⁻² and ABG, t ha⁻¹. CV stands for the coefficient of variation.

	Growth Stage	Min	Max	Mean	CV
LNC	Full Data set	1.41	5.16	2.61	38.63
	Tillering	4.00	5.16	4.65	7.69
	Jointing	2.43	3.30	2.84	8.93
	Booting	1.41	2.33	1.83	12.95
	Heading	1.77	2.32	2.07	6.52
PNC	Full Data set	0.83	4.15	1.71	56.29
	Tillering	3.09	4.15	3.68	8.91
	Jointing	1.70	2.37	2.00	9.87
	Booting	0.88	1.49	1.19	11.54
	Heading	0.83	1.20	1.01	9.90
LNA	Full Data set	1.12	6.41	3.48	33.04
	Tillering	1.26	2.93	2.05	24.21
	Jointing	2.89	6.41	4.29	22.48
	Booting	2.41	5.78	3.92	21.61
	Heading	1.12	4.23	2.64	27.31
PNA	Full Data set	18.21	156.36	77.34	39.58
	Tillering	18.21	39.00	28.01	21.83
	Jointing	46.75	95.54	65.96	19.91
	Booting	52.50	115.94	79.98	17.79
	Heading	55.61	141.83	95.72	20.64
LCC	Full Data set	22.38	39.88	30.63	12.44
	Tillering	33.78	39.88	37.36	4.55
	Jointing	28.67	35.46	31.57	5.09
	Booting	25.20	34.92	28.55	7.38
	Heading	27.69	34.37	30.26	5.08
CCC	Full Data set	29.36	293.72	118.45	51.89
	Tillering	32.71	81.63	56.20	23.99
	Jointing	65.51	201.15	122.81	27.85
	Booting	102.74	293.72	186.85	23.35
	Heading	29.36	127.69	86.03	26.14
ABG	Full Data set	0.54	12.37	6.10	54.24
	Tillering	0.54	0.98	0.75	14.86
	Jointing	2.56	4.19	3.28	13.54
	Booting	5.55	8.56	6.69	9.85
	Heading	6.69	11.99	9.36	13.43

depicts the variation in the measured LNC, PNC, LNA, PNA, LCC, CCC and ABG. In the full data set (with all growth stages combined), the LNC varies from 1.41% to 5.16% with coefficient of variation (CV) at 38.63%, the PNC varies from 0.83% to 4.15% with CV at 56.29%, the LNA varies from 1.12 g m⁻² to 6.41 g m⁻² with CV at 33.04%, the PNA varies from 18.21 kg ha⁻¹ to 156.36 kg ha⁻¹ with CV at 39.85%, the LCC varies from 22.38 μg cm⁻² to 39.88 μg cm⁻² with CV at 12.44%, the CCC ranges from 29.39 μg cm⁻² to 293.72 μg cm⁻² with CV at 51.98%, and the ABG varies from 0.54 t ha⁻¹ to 12.38 t ha⁻¹ with CV at 54.24%.

The means of LNC and PNC generally show a decreasing pattern from the tilling to heading growth stages, due to dilution effects [4]. Apart from that, the mean LNC at the heading growth stage is slightly higher than the mean at the booting growth stage. This is because our LNC is measured on the green leaves. During the heading growth stage, nitrogen is allocated from senescent leaves to young (green) leaves, which yields a LNC higher than that measured by combining senescent and green leaves together. The means of LCC show a pattern similar to those of LNC. The means of LNA and CCC first increase and then decrease from the tillering to heading growth stages. The means of PNA and ABG show an increasing pattern from the tillering to heading growth stages. The CV of each trait is lower for the single growth stage (data collection campaign) than for the full data set. The maximum values of the ratio between the single growth stage CV and full data set CV are 0.34, 0.20, 0.83, and 0.55 for LNC, PNC, LNA, and PNA, respectively.

Figure 3 shows the pairwise scatter plots and correlation coefficients of all seven traits (LNC, PNC, LNA, PNA, LCC, CCC and ABG). LNC is more correlated with LCC than CCC and ABG regardless of whether for the full dataset or a specific growth stage. PNC is highly correlated with LNC regardless of whether it is for the full data set or a specific growth stage. However, PNC is more correlated with LCC than CCC and ABG on the full data set, similar to LNC. However, at a specific growth stage, PNC is more correlated with CCC. Both LNA and PNA are more correlated with ABG than LCC and CCC for almost all specific growth stages. For the full data set, LNA is more correlated with CCC, and PNA is more correlated with ABG of the LCC, CCC and ABG traits. Except for LNC vs. PNC, all other pairwise relationships between LNC, PNC, LNA, PNA, LCC, CCC and ABG are affected by the growth stages. Although the growth-stage-specific correlation coefficients are positive, the full data set correlation coefficients of LNC vs. CCC, LNC vs. ABG, PNC vs. CCC, PNC vs. ABG, LNA vs. LCC, and PNA vs. LCC are negative.

3.2. Performance of the VIs

Figure 4 depicts the full-growth-stage and growth-stage-specific MRE_cv and R${}_{cv}^2$ of selected VIs for LNC, PNC, LNA and PNA prediction. For growth-stage-specific models of each tillering, jointing and heading growth stage, almost all selected VIs show good performance with relatively low MRE_cv and high R${}_{cv}^2$ values for all four nitrogen traits, except for TCARI and NDVI_483,503, whose performance is slightly worse than that of the other VIs. For the booting growth stage, the R${}_{cv}^2$ values of the growth-stage-specific models are lower than those of the other three growth stages for all four nitrogen traits. For the full-growth-stage models, the performance of all selected VIs is lower than the corresponding growth-stage-specific models (except for the booting growth stage).

The performance of the selected VIs is more diverse in the full-growth-stage models and booting growth stage growth-stage-specific models than in tillering, jointing and heading growth stage growth-stage-specific models. The MRE_cv values of the growth-stage-specific models for LNC and PNC are generally lower than those of the corresponding models for LNA and PNA. The MRE_cv values of the full-growth-stage models for LNA and PNA are generally lower than those of the corresponding models for LNC and PNC.

Table 3. The goodness-of-fit of the best-performing VIs of the full-growth-stage and growth-stage-specific models for LNC, PNC, LNA and PNA selected according to MRE_cv.

			Calibration			Validation
	Data Set	Vegetation Index	MREcv	RMSEcv	R${}_{\mathbf{cv}}^2$	MREval	RMSEval	R${}_{\mathbf{val}}^2$
LNC	Full Data set	TCARI/OSAVI	25.85	0.74	0.62	24.95	0.68	0.61
	Tillering	CIgreen	3.34	0.17	1.00	2.03	0.11	0.88
	Jointing	REPLi	3.80	0.12	1.00	4.06	0.13	0.87
	Booting	CIre	10.00	0.20	0.50	10.48	0.21	0.08
	Heading	TCARI/OSAVI	3.87	0.09	0.94	4.31	0.14	0.07
PNC	Full Data set	TCARI/OSAVI	30.85	0.62	0.69	25.18	0.53	0.73
	Tillering	CIgreen	3.45	0.14	1.00	2.25	0.10	0.90
	Jointing	mND705	3.20	0.08	1.00	2.50	0.06	0.95
	Booting	NDPI	6.06	0.09	0.65	5.09	0.08	0.50
	Heading	mND705	4.38	0.05	1.00	6.86	0.08	0.39
LNA	Full Data set	TBVITian	16.91	0.61	0.72	17.82	0.74	0.69
	Tillering	D735	8.13	0.18	1.00	7.05	0.20	0.82
	Jointing	OSAVI2	11.76	0.56	0.99	9.78	0.47	0.81
	Booting	NDPI	14.10	0.63	0.61	14.49	0.63	0.21
	Heading	NDVI	16.24	0.52	0.99	15.01	0.45	0.28
PNA	Full Data set	TCARI/OSAVI	12.43	10.69	0.84	16.75	12.13	0.80
	Tillering	D735	7.52	2.26	1.00	6.78	2.27	0.85
	Jointing	NDRE	10.38	7.82	1.00	9.98	6.70	0.79
	Booting	TCARI2/OSAVI2	8.19	7.75	0.74	9.53	10.75	0.36
	Heading	REPLE	10.30	11.44	1.00	13.35	13.94	0.51

lists the best-performing VI according to the MRE_cv values for each full-growth-stage and growth-stage-specific model of each nitrogen trait. Figure 5 shows the scatter plot between the observed and predicted nitrogen traits in the corresponding validation data set of these models. For full-growth-stage models, the transformed chlorophyll absorption in the reflectance index/optimized soil-adjusted vegetation index (TCARI/OSAVI) yields the lowest MRE_cv values for LNC and PNC and PNA, while the three band vegetation index (TBVI_Tian) yields the lowest MRE_cv for LNA. The validation MRE (MRE_val) values are 24.95%, 25.18%, 17.82% and 16.75% and the validation R² (R${}_{val}^2$) values are 0.61, 0.73, 0.69 and 0.80 for LNC, PNC, LNA and PNA, respectively. The full-growth-stage models of LNA and PNA could yield lower MRE and higher R² values than LNC and PNC both in model calibration and validation.

For growth-stage-specific models, the MRE and R² range from 3.34% to 10.48% and 0.08 to 1.00 for LNC in both model calibration and validation. According to MRE_cv, the model performances are ordered as follows: tillering > jointing > heading > booting. The MRE and R² range from 3.20% to 6.86% and 0.39 to 1.00 for PNC in both model calibration and validation. The model performance order is similar to that of LNC. The validation MRE and R² range from 7.05% to 16.24% and 0.21 to 1.00 for LNA in both model calibration and validation. According to MRE_cv, the model performances are ordered as follows: tillering > jointing > booting > heading. The MRE and R² range from 6.78% to 13.35% for PNA in both model calibration and validation. According to MRE_cv, the model performances are ordered as follows: tillering > booting > heading > jointing. For each nitrogen trait, the growth-stage-specific model always yields a lower MRE in both model calibration and validation than the corresponding full-growth-stage models.

3.3. Performance of the PLSR and ML Methods

Figure 6 shows the growth-stage-specific and full-growth-stage MRE_cv and R${}_{cv}^2$ values of the PLSR and three ML methods (the ANN, RF and SVM) for the four nitrogen traits (LNC, PNC, LNA and PNA). There are no obvious patterns showing a certain method clearly outperforming the other methods. However, the PLSR may yield a large MRE_cv in the full-growth-stage models. Of the four nitrogen traits, the booting growth stage models yield the lowest R_cv and moderate or highest MRE_cv of the growth-stage-specific models and full-growth-stage models.

For LNC and PNC, the full-growth-stage models yield R${}_{cv}^2$ values comparable to those of the tillering, jointing and heading growth-stage-specific models. For LNA and PNA, the full-growth-stage models yield slightly lower R${}_{cv}^2$ values than the tillering, jointing and heading growth-stage-specific models. However, the R${}_{cv}^2$ values are still above 0.60. The MRE_cv values of LNC and PNC are slightly lower than the corresponding MRE values of LNA and PNA.

Table 4. The goodness-of-fit values of the best-performing models of the full-growth-stage and each growth-stage-specific models for LNC, PNC, LNA and PNA selected by MRE_cv.

			Calibration (Cross-Validated)			Validation
	Data Set	Method	MREcv	RMSEcv	${\mathit{R}}_{\mathbf{cv}}^2$	MREval	RMSEval	${\mathit{R}}_{\mathbf{val}}^2$
LNC	Full Data set	ANN	7.82	0.21	0.95	6.62	0.19	0.97
	Tillering	RF	3.10	0.16	0.98	3.38	0.19	0.86
	Jointing	SVM	3.86	0.12	1.00	6.80	0.20	0.44
	Booting	PLSR	11.42	0.23	0.15	8.76	0.19	0.50
	Heading	SVM	3.66	0.09	0.94	10.54	0.25	0.01
PNC	Full Data set	ANN	6.41	0.12	0.98	5.15	0.11	0.99
	Tillering	PLSR	3.02	0.12	0.98	2.56	0.12	0.90
	Jointing	SVM	3.03	0.07	1.00	4.65	0.11	0.80
	Booting	PLSR	7.30	0.10	0.44	3.84	0.06	0.82
	Heading	PLSR	5.16	0.06	0.95	6.59	0.07	0.83
LNA	Full Data set	SVM	14.85	0.60	0.72	11.04	0.49	0.81
	Tillering	ANN	7.10	0.15	0.98	9.18	0.20	0.91
	Jointing	PLSR	11.09	0.53	0.98	10.88	0.45	0.84
	Booting	ANN	14.59	0.67	0.54	15.73	0.83	0.22
	Heading	PLSR	14.65	0.48	0.99	28.26	0.55	0.88
PNA	Full Data set	SVM	11.32	9.16	0.84	11.11	11.40	0.87
	Tillering	PLSR	6.02	1.79	1.00	9.31	2.80	0.95
	Jointing	RF	9.79	6.88	0.99	9.32	7.54	0.69
	Booting	PLSR	8.59	8.85	0.61	6.10	7.50	0.77
	Heading	PLSR	10.71	11.66	0.90	10.83	9.49	0.95

lists the best-performing models of the full-growth-stage and each growth-stage-specific models for LNC, PNC, LNA and PNA selected by MRE_cv. Figure 7 shows the scatter plots between the observed and predicted nitrogen traits in the corresponding validation data sets. In the full-growth-stage models, the ANN method outperforms the other three methods for LNC and PNC, while the SVM outperforms the other three methods for LNA and PNA. For both LNC and PNC, the full-growth-stage models built on the ANN yield MRE below 7.82% and R² above 0.95 both in model calibration and validation. For LNA and PNA, the performance of the full-growth-stage models slightly degrades. The MRE ranges between 11.04% and 14.85%, and R² ranges between 0.72 and 0.87 in both model calibration and validation.

For LNC, only the tillering growth-stage-specific model yields good performance, with MRE below 3.38% and R² above 0.86 in both model calibration and validation. Although the growth-stage-specific models of the other three growth stages (jointing, booting and heading) yield lower MRE_val, the R${}_{val}^2$ values are lower (ranging between 0.01 and 0.50). For LNA, all tillering, jointing and heading growth-stage-specific models yield good performance, with MRE below 28.26% and R² above 0.84 in both model calibration and validation. The booting growth-stage-specific model yields lower R${}_{val}^2$ (R${}_{val}^2$ = 0.22). For both PNC and PNA, the growth-stage-specific models from tillering to heading growth stages all show good performance, with MRE below 10.83% and R² above 0.44 in both model calibration and validation.

Comparing Table 3 and Table 4, which present the best-performing VI-based and full-spectra-based models, respectively, the latter is clearly superior, regarding the full-growth-stage models of each nitrogen trait. The VI-based growth-stage-specific model of the booting growth stage yields lower accuracy for each nitrogen trait. With the full-spectra-based models, the prediction accuracy is improved for PNC and PNA.

3.4. Evaluation on Hyperspectral Imageries

Figure 8 shows the pixelwise inversion results of the four N traits (LNC, PNC, LNA and PNC), using the corresponding optimal ML full-growth-stage models listed in Table 4 (the ANN for LNC and PNC and the SVM for LNA and PNA). The N treatment level variation is relatively homogenous. The spatial pattern due to different nitrogen treatments is more obvious in the maps of LNA and PNA than those of LNC and PNC. The N trait variations caused by growth stages are in line with the observations (detailed in Table 2).

4. Discussion

With a focus on fresh leaves or plant canopies instead of dry leaves or dried and ground leaves, the remote sensing-based estimation of nitrogen traits is mainly based on the relationship between chlorophyll and nitrogen by means of the use of visible and near infrared (VNIR) spectra. The pairwise correlation analysis (Figure 3) shows that the relationships among the four nitrogen traits and those between the four nitrogen traits and the three biochemical traits (LCC, CCC and ABG) are obviously affected by the growth stage, except for the relationship between LNC and PNC. During the vegetative growth stage, N is abundant in the vegetative tissues, while in the reproductive growth stage, N reallocates from the leaves to the reproductive organs (ears, fruits or seeds). In the reallocation process, the assembly and disassembly of N are accompanied by canopy structure variation. Additionally, the ABG increases more slowly with growth stage development, and the “dilution effects” end at the heading growth stage when the plant N content dominates the canopy reflectance [64]. The reallocation of N and variation in the biomass increase the rate change of the relationships between nitrogen and chlorophyll or ABG across the growth stages.

The performance of the growth-stage-specific VI-based nitrogen trait estimation models is relatively constant within the same growth stage. A similar pattern was observed in a previous study [15]. Patel et al. [15] evaluated the performance of several VIs for irrigated perennial ryegrass PNC estimation and found that the selected VIs showed a similar correlation with PNC when applied to single growth stage data. When multiple growth stage data were pooled together, the accuracy of VI-based full-growth-stage models degraded obviously. Li et al. [16] and Patel et al. [15] found similar results. This is because, in addition to the nitrogen and chlorophyll content, canopy reflectance is affected by LAI, ABG, nonphotosynthetic components, soil background and other factors. With the growth stage development, the canopy structure and background situation vary. The VIs that use several bands cannot fully capture the intrinsic relationship between the canopy spectrum and the target variable across growth stages.

There is general interest in exploring the red-edge spectrum, which are bands where reflectance sharply increases from the lower reflectance red region to high reflectance near the infrared region, for biochemical parameter estimation [13,25,28,64,65]. The REP, which is the main inflection point of the red-edge spectrum, has been found to be useful for chlorophyll and nitrogen content estimation. In this work, the REPs calculated with different techniques, such as linear four-point interpolation, linear extrapolation and Gaussian fit, were evaluated. Despite taking advantage of the narrower continuing hyperspectral data, the REPs are not obviously superior to other VIs in either growth-stage-specific or full-growth-stage models. However, there are several VIs that use the red-edge band—such as CI_re, mND705, NDPI, and D735—that perform better than other VIs. Several works have pointed out the importance of blue bands for dry-mass-based nitrogen content (LNC and PNC) estimation [16,19,24]; however, our sensor does not cover this region.

Compared to the VI-based models, the PLSR and ML models generally yield better performance for the full-growth-stage models. This may be because more bands are involved, and the full potential of the hyperspectral data is explored. For the growth-stage-specific models, the PLSR and ML models generally show similar performance to, or even worse performance than, the corresponding VI-based models. This may indicate that for growth-stage-specific models, the VI-based estimation is optimal and cannot be further improved.

5. Conclusions

This work evaluates the performance of univariate regression model vegetation indices (VIs), traditional multivariate calibration method partial least squares regression (PLSR) and modern machine learning (ML) methods, such as the artificial neural network (ANN), random forest (RF), and support vector machine (SVM), for paddy rice nitrogen trait—leaf nitrogen content (LNC), plant nitrogen content (PNC), leaf nitrogen accumulation (LNA) and plant nitrogen accumulation (PNA)—estimation, with a special focus on the effects of growth stages. All methods were evaluated both throughout the growing season and at each growth stage (tillering, jointing, booting and heading).

The results indicate that the correlations between the four nitrogen traits and the other three biochemical traits—leaf chlorophyll content (LCC), canopy chlorophyll content (CCC) and aboveground biomass (ABG)—are affected by the growth stage. Within a single growth stage, the performance of the selected VIs is relatively constant. For full-growth-stage models, the performance of the VI-based models is more diverse. For full-growth-stage models, the TCARI/OSAVI performed best for LNC, PNC and PNA estimation, while the TBVI_Tian performed best for LNA estimation. There are no obvious patterns showing whether the PLSR, ANN, RF or SVM performs best in either the growth-stage-specific or full-growth-stage models.

For the growth-stage-specific models, a lower MRE and higher R² can be acquired at the tillering and jointing growth stages. The PLSR and ML methods yield obviously better estimation accuracy for the full-growth-stage models than the VI-based models. For the growth-stage-specific models, the performance of VI-based models seems optimal and cannot be obviously improved. These results suggest that building linear regression models on VIs for paddy rice nitrogen traits estimation is still a reasonable choice when only a single growth stage is involved. However, when multiple growth stages are involved or the phenology information is missing, using PLSR or ML methods is a better option. Nevertheless, it must be noticed that only one paddy rice variety was tested in our study. Whether our conclusions still hold when other varieties/cultivations or multiple varieties/cultivations are involved needs further evaluation.