Enhancing the Nitrogen Signals of Rice Canopies across Critical Growth Stages through the Integration of Textural and Spectral Information from Unmanned Aerial Vehicle (UAV) Multispectral Imagery

Abstract

This paper evaluates the potential of integrating textural and spectral information from unmanned aerial vehicle (UAV)-based multispectral imagery for improving the quantification of nitrogen (N) status in rice crops. Vegetation indices (VIs), normalized difference texture indices (NDTIs), and their combination were used to estimate four N nutrition parameters leaf nitrogen concentration (LNC), leaf nitrogen accumulation (LNA), plant nitrogen concentration (PNC), and plant nitrogen accumulation (PNA). Results demonstrated that the normalized difference red-edge index (NDRE) performed best in estimating the N nutrition parameters among all the VI candidates. The optimal texture indices had comparable performance in N nutrition parameters estimation as compared to NDRE. Significant improvement for all N nutrition parameters could be obtained by integrating VIs with NDTIs using multiple linear regression. While tested across years and growth stages, the multivariate models also exhibited satisfactory estimation accuracy. For texture analysis, texture metrics calculated in the direction D3 (perpendicular to the row orientation) are recommended for monitoring row-planted crops. These findings indicate that the addition of textural information derived from UAV multispectral imagery could reduce the effects of background materials and saturation and enhance the N signals of rice canopies for the entire season.

1. Introduction

As one of the most important staple crops around the world, rice (Oryza sativa L.) crop feeds more than 50% of the world’s population. Nitrogen is the primary element for crop growth due to the crucial impact on crop yield formulation and grain quality determination. Monitoring N status in rice leaves or plants could provide valuable information for growth diagnosis [1] and precise field N management [2], so as to improve the N use efficiency and reduce environmental pollution.

Crop N status is often represented by LNC and PNC in precision agriculture. However, the majority of studies reported that it was difficult to quantify LNC and PNC precisely with the canopy reflectance data and N status monitoring was strongly influenced by growth stage [3,4]. Nitrogen signals from rice leaves and plants captured by canopy spectra at early growth stages are weak, because the background materials (e.g., soil, water) occupy a large proportion of the field of view and the biomass of rice crops increases as the N concentration decreases due to the N dilution effect [5,6]. Since those negative influences could be reduced at reproductive stages with the stabilization of leaf biomass [6], the N status is difficult to detect at early stages and could be estimated more accurately for late stages. Zhou et al. [6] found that the LNC at vegetative stages (with low canopy coverage) should be estimated with stage-specific regression models, and only the LNC at reproductive stages (with high canopy coverage and even canopy closure) could be estimated with a common model. In contrast to the influence of background materials for early stages, the low sensitivity of spectral signals to the growth dynamics of dense canopies becomes another issue for the estimation of N concentration for late stages (from booting to grain filling). Therefore, determining how to enhance N signals and build suitable models for N status monitoring over the whole season remains to be addressed.

Other than spectral information, the textural information inherent in the UAV imagery has recently been proved to have great potential in crop growth monitoring [7,8]. Image texture analysis can be applied to the identification of image tonal variation through a function of local variance in the images [9]. The texture metrics derived from spectral bands can characterize canopy structural patterns and enhance the sensitivity of remotely sensed data to biophysical properties, especially at dense canopy coverage conditions [10]. During the early growth period (from tillering to booting stages) of rice crops, the leaves often exhibit remarkable variation in color as a result of N nutrition consumption for vegetative development and the topdressing of N fertilizer [11,12]. These tonal variations within and between crop canopies could be captured by texture metrics. Therefore, the image-derived texture metrics have the potential to better characterize the N status variation at critical growth stages. The availability of different data sources enables data fusion, i.e., the combination of datasets from two or more data types with different characteristics [7]. Data fusion may allow a more holistic interpretation of the relationship between remote sensing information and crop parameters. For example, the fusion of spectral information and textural information from UAV imagery showed an increased prediction performance for crop biomass [8]. Since the textural and spectral features respond to the N status dynamics in different ways, the integration of their complementary information could be useful for improving the estimation of N nutrition parameter across critical growth stages. The performance of texture analysis depends on the choices of texture metrics, the spectral band used for texture derivation, window size, and texture direction [13,14]. While the former three factors were well studied for vegetation monitoring, the effect of texture direction has rarely been investigated. Previous studies often conducted texture analysis with the default direction (45°), because most of them focused on naturally grown forests with disorderly distributed tress [7,15]. Others exercised with different directions for naturally grown trees, but they did not explicitly provide the mechanism underlying the direction choices [10,14,16]. Rice crops are mostly planted in rows and exhibit regular spatial patterns, but how the texture direction affects the performance of derived textural information remains unclear. Whether the performance of across-row texture metrics differs from that of along-row texture metrics and whether those two directions outperform the default diagonal direction need to be investigated.

Therefore, our aims were (1) to determine the optimal texture metrics and texture indices from UAV multispectral imagery for estimating rice N nutrition parameters, (2) to investigate the directional effect of texture analysis on N status estimation for row-planted rice crops, and (3) to examine the feasibility of integrating spectral and textural information in the improvement of N status monitoring.

2. Materials and Methods

2.1. Experimental Design

In 2015 and 2016, we conducted two field experiments in Rugao, Jiangsu Province, China (120°45′E, 32°16′N) [8]. In 2015, there were four N rates (0, 100, 200, and 300 kg N ha⁻¹). The minimum and maximum N rates had at one density (22 × 10⁴ plants ha⁻¹) and the intermediate N rates had two densities (13 × 10⁴ plants ha⁻¹ and 22 × 10⁴ plants ha⁻¹). The ratio of N fertilizer application was 4:1:3:2 for basal fertilizer before transplanting, tillering stage, jointing stage, and booting stage. In 2016, there were three N rates (0, 150, and 300 kg N ha⁻¹) and two densities (13 × 10⁴ plants ha⁻¹ and 22 × 10⁴ plants ha⁻¹). The ratio of N fertilizer application was 4:2:2:2 for basal fertilizer before transplanting, tillering stage, jointing stage, and booting stage. There were two rice cultivars, (Wuyunjing 24, V1) and (Y Liangyou 1, V2). Each experiment had 36 plots (12 treatments with three replications) with 30 m² (6 m × 5 m).

2.2. Ground Sampling and UAV Image Acquisition

The UAV campaigns and ground destructive samplings were taken at critical growth stages of rice (

Table 1. Experimental designs and dates of data collection on the ground.

Year	2015	2016
N rate (kg ha−1)	0 (N0), 100 (N1), 200 (N2), 300 (N3)	0 (N0), 150 (N1), 300 (N2)
Planting density (plants ha−1)	13 × 104, 22 × 104	13 × 104, 22 × 104
UAV flight date	5 August, 14 August, 26 August, 9 September	6 August, 14 August, 28 August, 9 September
Sampling date	31 July, 15 August, 26 August, 10 September	6 August, 14 August, 28 August, 8 September
Growth stage	Jointing, Early-booting, Late-booting, Grain filling	Jointing, Booting, Heading, Grain filling

). After the UAV campaign, three hills of rice plants were randomly harvested in each plot and separated into different organs. All the samples were oven dried at 105 °C for 30 mins, followed by drying at 80 °C, until a constant weight was reached. Then, they were weighed, ground, and stored in plastic bags for laboratory chemical analysis. Aboveground biomass (AGB) was the total dry weight of all plant organs per unit ground area. Total N concentrations (%) were determined with the micro-Keldjahl analysis. The LNC and PNC were expressed on a dry weight basis (%), and LNA and PNA were expressed as N mass per unit ground area (g m⁻²). The LNA (g m⁻²) was obtained through multiplying LNC (%) by leaf dry biomass (LDB, t ha⁻¹). The PNA (g m²) was obtained through multiplying PNC (%) by AGB (t ha⁻¹).

In this study, we used an eight-rotor aircraft (Mikrokopter OktoXL) which had a maximum payload capacity of 2.5 kg and a flight duration of 10−20 min, depending on the battery and actual payload. A multispectral camera (Mini-MCA6; Tetracam, Inc., Chatsworth, CA, USA) was mounted on the UAV to acquire images in this study (Table 1). The specification of this multispectral camera can be found in Reference [17]. The UAV was flown at 100 m above ground level, resulting in a nominal resolution of 0.054 m in the multispectral images. The camera was set to continuous data capture at one frame per three seconds with a fixed aperture and exposure according to the light conditions. Aerial images were saved in the camera memory card in a 10 bit RAW format.

The UAV campaigns were conducted under a clear sky and calm wind conditions between 11:00 and 13:30 (local time). After each flight, only one image that covered all 36 plots in the nadir position was selected for subsequent analysis.

2.3. Image Processing

The UAV-based image preprocessing workflows were performed by following Reference [17] and processed in IDL/ENVI environment (Exelis Visual Information Solutions, Boulder, Colorado, USA). The pre-processing workflows included noise reduction, vignetting correction, and lens distortion correction. Subsequently, the six bands were co-registered with 25 ground control points (GCPs) evenly distributed in the study area (Figure 1). Then, we conducted radiometric correction with an empirical line correction method using a six flat calibration canvas at different reflectance intensities [17,18].

The reflectance value of each plot was extracted as the mean over a region of interest (ROI) in the non-sampling area. Then, ten VIs commonly used in crop growth monitoring were calculated from the reflectance values (

Table 2. List of vegetation indices (VIs) used in this study.

	Index	Formula	Reference
Color indices	Normalized green red difference index (NGRDI)	$NGRDI=(R_{550}-R_{680})/(R_{550}+R_{680})$	[20]
	Visible atmospherically resistance index (VARI)	$VARI=(R_{550}-R_{680})/(R_{550}+R_{680}-R_{490})$	[21]
	Green leaf index (GLI)	$GLI=(2R_{550}-R_{680}-R_{490})/(2R_{550}+R_{680}+R_{490})$	[22]
RE indices	MERIS terrestrial chlorophyll index	$MTCI=(R_{800}-R_{720})/(R_{720}+R_{680})$	[23]
	red edge chlorophyll index (CIRE)	$C{I}_{RE}=R_{800}/R_{720}-1$	[24]
	Normalized difference red edge index (NDRE)	$NDRE=(R_{800}-R_{720})/(R_{800}+R_{720})$	[25]
NIR indices	Normalized difference vegetation index (NDVI)	$NDVI=(R_{800}-R_{680})/(R_{800}+R_{680})$	[26]
	Renormalized difference vegetation index (RDVI)	$RDVI=(R_{800}-R_{680})/\sqrt{(R_{800}+R_{680})}$	[27]
	Optimized soil adjusted vegetation index (OSAVI)	$OSAVI=(1+0.16)(R_{800}-R_{680})/(R_{800}+R_{680}+0.16)$	[28]
	Modified triangular vegetation index 2 (MTVI2)	$MTVI2=\frac{1.5[1.2(R_{800}-R_{550})-2.5(R_{670}-R_{550})]}{\sqrt{{\left(2R_{800}+1\right)}^2-\left(6R_{800}-5\sqrt{R_{670}}\right)}-0.5}$	[29]

The bands used in the formula corresponded to the bands in the multispectral imagery.

).

The grey level co-occurrence matrix (GLCM) was employed in this study to evaluate the potential of texture analysis on reflectance images for improving rice N status monitoring. Eight GLCM-based texture metrics including mean (MEA), variance (VAR), homogeneity (HOM), contrast (CON), dissimilarity (DIS), entropy (ENT), second moment (SEM), and correlation (COR) were computed using the ENVI software. A detailed description of the eight texture metrics can be found in Reference [19]. Since the smallest window size (3 × 3) was comparable to the row spacing at jointing stage and larger window sizes did not exhibit significant differences (data not shown), only the window size (3 × 3) was investigated. Next, inter-correlations between different texture metrics were analyzed to decrease the data dimensionality and improve the data processing efficiency. Finally, we chose four texture metrics (i.e., MEA, CON, DIS, and COR) at the minimal window size (3 × 3) for texture analysis with five bands (blue, green, red, RE, and NIR) in four directions (D1 = 0°, D2 = 45°, D3 = 90° and D4 = 135°).

Based on individual texture bands, the normalized difference texture index (NDTI) was also derived to improve the performance of texture analysis. The NDTI was calculated in the equation below:

$$\mathrm{NDTI}(T_1,T_2)=\frac{(T_1-T_2)}{(T_1+T_2)}$$

(1)

where T₁ and T₂ are a random texture feature from the five bands in four directions.

2.4. Statistical Analysis

The data collected from the two field experiments were pooled to examine the relationships of agronomic variables with VIs, NDTIs, and their combinations by simple regression (SR) and stepwise multiple linear regression (SMLR). Model validation was performed on the pooled data using a k-fold (k = 10) cross-validation procedure and evaluated by the root mean square error (RMSE) and the relative RMSE (RRMSE).

3. Results

3.1. Relationships of N Nutrition Parameters with VIs

Table 3. Relationships (R²) among different N nutrition parameters.

	LDB	AGB	LNC	LNA	PNC	PNA
LDB	1.00
AGB	0.55	1.00
LNC	0.20	0.00	1.00
LNA	0.91	0.35	0.46	1.00
PNC	0.13	0.02	0.78	0.34	1.00
PNA	0.81	0.71	0.19	0.77	0.13	1.00

LDB is leaf dry biomass, AGB is aboveground biomass, LNC is leaf nitrogen concentration, LNA is leaf nitrogen accumulation, PNC is plant nitrogen concentration, PNA is plant nitrogen accumulation. Bold numbers indicate values over 0.70.

shows the relationships among different agronomic parameters. The LNC and PNC and LNA and PNA were highly correlated to each other. In addition, biomass was highly correlated to both LNA and PNA.

Table 4. Coefficient of determination (R²) values for the relationships between VIs and N nutrition parameters in rice.

	LNC (%)		LNA (g m−2)		PNC (%)		PNA (g m−2)
VI	L	E	L	E	L	E	L	E
NGRDI	0.00ns	0.00ns	0.24***	0.27***	0.00ns	0.01ns	0.33***	0.40***
VARI	0.00ns	0.00ns	0.32***	0.34***	0.00ns	0.00ns	0.33***	0.41***
GLI	0.01ns	0.01ns	0.12***	0.15***	0.03*	0.05***	0.26***	0.32***
MTCI	0.17***	0.16***	0.73***	0.68***	0.12***	0.10***	0.66***	0.62***
CIRE	0.20***	0.20***	0.77***	0.72***	0.15***	0.13***	0.68***	0.65***
NDRE	0.20***	0.19***	0.74***	0.77***	0.14***	0.13***	0.68***	0.73***
NDVI	0.03*	0.03*	0.39***	0.47***	0.02*	0.01ns	0.44***	0.63***
RDVI	0.07***	0.07***	0.44***	0.50***	0.06***	0.05***	0.60***	0.72***
OSAVI	0.07***	0.06***	0.45***	0.52***	0.06***	0.05***	0.58***	0.73***
MTVI2	0.06***	0.06***	0.45***	0.46***	0.05***	0.04**	0.59***	0.66***

The numbers in bold correspond to the best-fit function for each N nutrition parameter. Significance level: “ns” represents no significance, * p < 0.01, ** p < 0.001, *** p < 0.0001. “L” and “E” indicate linear and exponential regression.

shows the relationships between N nutrition parameters and VIs derived from UAV multispectral images. For foliar parameters, the relationships between VIs and LNC were remarkably low, and the highest R² was only 0.2 obtained by CI_RE and NDRE (Figure 1a). These two indices also exhibited the strongest relationships with LNA (R² = 0.77 for both) (Figure 1b). For plant parameters, PNC was most correlated to CI_RE as seen for LNC but with a lower R² value (R² = 0.15) (Figure 1c). The OSAVI exhibited the strongest relationships with PNA (R² = 0.73) which was also found for NDRE (R² = 0.73) (Figure 1d). In general, RE-based VIs exhibited better relationships with N nutrition parameters than other indices.

3.2. Texture Analysis and Performance of Texture Indices

3.2.1. Relationships between N Nutrition Parameters and Individual Texture Metrics

Figure 2 shows the relationships between N nutrition parameters and single-band texture metrics. For LNC and PNC estimation, the texture metrics in the direction D3 were generally superior to those in other directions, while directions did not influence the performance of the optimal texture feature for LNA and PNA estimation. The DIS feature and the MEA feature exhibited a higher correlation than other texture metrics with N concentration and N accumulation indicators, respectively. The DIS feature in the RE band was superior to other bands in LNC estimation with an R² of only 0.26 (Figure 2a). The MEA feature in the red band was most related to LNA (R² = 0.39) among all the texture metrics (Figure 2b). Similar to LNC, the DIS feature in the RE band exhibited a higher correlation with PNC compared to their counterparts (Figure 2c). The MEA feature in the NIR band obtained a higher correlation with PNA (R² = 0.55) than other texture metrics (Figure 2d).

3.2.2. Relationships between N Nutrition Parameters and Texture Indices

Considerable improvement was obtained with texture indices for N nutrition parameters estimation comparing with individual texture metrics (

Table 5. Relationships (R²) between N nutrition parameters and the top four best-performing texture indices.

	NDTI1	NDTI2	NDTI3	NDTI4
LNC (%)	(DISB5_D4, CONB4_D3)	(DISB4_D2, CORB5_D3)	(DISB4_D3, CORB5_D2)	(DISB4_D3, CORB5_D3)
R2	0.31	0.30	0.30	0.29
LNA (g m−2)	(MEAB5_D3, MEAB4_D4)	(MEAB5_D3, MEAB4_D3)	(MEAB5_D4, MEAB4_D4)	(MEAB5_D4, MEAB4_D3)
R2	0.72	0.72	0.72	0.72
PNC (%)	(DISB4_D3, CORB5_D3)	(DISB5_D4, CONB4_D3)	(CONB4_D3, CORB5_D2)	(DISB4_D4, CORB5_D3)
R2	0.41	0.38	0.37	0.36
PNA (g m−2)	(MEAB5_D3, MEAB2_D4)	(MEAB5_D4, MEAB2_D4)	(MEAB5_D3, MEAB2_D3)	(MEAB5_D4, MEAB2_D3)
R2	0.75	0.75	0.75	0.75

All regressions are statistically significant (p < 0.0001).

). The majority of optimal NDTIs were composed of texture metrics in D3 and D4. The optimal NDTIs for LNC and PNC estimation consisted of DIS, CON, and COR features calculated from RE and NIR bands. While the optimal NDTIs for LNA estimation consisted of MEA features calculated from RE and NIR bands, and those for PNA estimation were composed of MEA features calculated from green and NIR bands. The index (NDTI(DIS_{B5_D4}, CON_{B4_D3})) was closely related to both LNC and PNC with R² of 0.29 and 0.41 (Figure 3a,c). The NDTI(MEA_{B5_D3}, MEA_{B4_D4}) performed best in the LNA estimation (R² = 0.72) with exponential regression (Figure 3b). The NDTI(MEA_{B5_D3}, MEA_{B2_D3}) performed best in the PNA estimation (R² = 0.75), and obvious saturation occurred at high PNA levels (Figure 3d).

3.3. Performance of the Integration of VIs and Texture Indices

Table 6. Stepwise multiple linear regression output for the relationship between each of the four N nutrition parameters and the best combination of vegetation and textural indices.

	R2	Variables and Intercept	B	p	Tol	VIF
LNC (%)	0.70	Intercept	1.022	0.00
		NDTI(DISB4_D2, CORB5_D3)	−0.450	0.00	0.931	1.074
		NDRE	2.256	0.00	0.160	6.264
		NDTI(MEAB5_D4, MEAB2_D4)	−1.392	0.00	0.161	6.217
LNA (g m−2)	0.85	Intercept	0.173	0.00
		NDRE	5.044	0.00	0.332	3.009
		NGRDI	−1.120	0.00	0.320	3.121
		NDTI(DISB4_D3, CORB5_D3)	−0.464	0.00	0.772	1.296
PNC (%)	0.68	Intercept	0.497	0.00
		NDTI(DISB4_D3, CORB5_D3)	−0.646	0.00	0.779	1.284
		NDTI(MEAB5_D3, MEAB4_D3)	2.316	0.00	0.210	4.757
		NDTI(MEAB5_D4, MEAB2_D4)	−1.485	0.00	0.217	4.604
PNA (g m−2)	0.86	Intercept	−0.498	0.00
		NDTI(MEAB5_D3, MEAB2_D4)	3.454	0.00	0.274	3.644
		MTVI2	2.621	0.00	0.405	2.470
		NGRDI	−1.850	0.00	0.296	3.384

VIF and Tol indicate variance inflation factor and tolerance.

shows the performance of the combination of VIs and texture indices in the estimation of N nutrition parameters with stepwise multiple linear regression. The hybrid model with VIs and NDTIs could explain 70% of the variability in LNC which was significantly higher than using VIs or texture indices alone. Compared with the optimal VI (NDRE) for LNA estimation, the combination of VIs and NDTIs produced a slightly better relationship. Similar to the LNC estimation, the multivariate model improved the regression significantly from 0.41 to 0.68 with three texture indices. Furthermore, the hybrid model explained more than 80% of the variability in PNA (R² = 0.86) which was higher than using individual remote sensing variables alone. Interestingly, all the multivariate models for plant-level variables included NDTI(MEA_{B5_D3}, MEA_{B2_D4}). Furthermore, all variables in the hybrid models were statistically significant and multicollinearity effects were not observed.

3.4. Model Validation

Table 7. Accuracy assessment with RMSE and RRMSE values for the estimation of N nutrition parameters with different predictor variables.

Predictor Variable	LNC (%)		LNA (g m−2)		PNC (%)		PNA (g m−2)
	RMSE	RRMSE	RMSE	RRMSE	RMSE	RRMSE	RMSE	RRMSE
VI	0.35	15.28%	1.35	24.56%	0.31	24.78%	3.16	32.20%
NDTI	0.32	14.24%	1.54	27.91%	0.26	20.67%	3.33	33.90%
NDTI+VI	0.22	9.70%	1.16	21.04%	0.19	15.34%	2.38	24.17%

shows the accuracy assessment for N nutrition parameters with the best-fit model in each category. The multivariate model yielded the significantly higher estimation accuracy (RMSE = 0.22 and RRMSE = 9.70%) for LNC estimation than using VIs and NDTIs alone (Figure 4a–c). Compared with the optimal VI, the combination of VI and NDTI exhibited a slightly higher estimation accuracy for LNA (RMSE = 1.16 g m⁻² and RRMSE = 21.04%) estimation. The highest estimation accuracy (RMSE = 0.19 and RRMSE = 15.34%) for PNC estimation was obtained by the hybrid model, which was significantly higher than other techniques (Figure 4d–f). Similarly, the multivariate models yielded the highest estimation accuracy for PNA (RMSE = 2.38 g m⁻² and RRMSE = 24.17%) estimation when compared with other techniques. Furthermore, the optimal multivariate models for LNC and PNC in Table 6 exhibited satisfactory performance when validated with different dataset groups (

Table 8. Accuracy assessment with RMSE and RE values for the optimal multivariate models of LNC and PNC on different dataset groups.

Year	Stage	LNC (%)		PNC (%)
		RMSE	RRMSE	RMSE	RRMSE
2015	Jointing	0.27	11.03%	0.27	19.95%
	Early-booting	0.22	9.48%	0.20	14.09%
	Late-booting	0.25	10.95%	0.19	15.65%
	Filling	0.18	9.45%	0.12	14.58%
	All	0.23	10.41%	0.20	16.82%
2016	Jointing	0.17	6.95%	0.18	11.55%
	Booting	0.26	10.83%	0.20	14.27%
	Heading	0.17	7.15%	0.22	17.80%
	Filling	0.22	10.72%	0.12	11.21%
	All	0.21	8.99%	0.18	13.93%

).

When all the texture indices in the multivariate models used the texture metrics calculated in the direction D3 for simplicity, the validation accuracies for the majority of N nutrition parameters were marginally lower (data not shown). For example, the validation RMSE and RRMSE values of the LNC estimation decreased slightly to 0.24% and 10.52%, respectively.

4. Discussion

4.1. Differences among N Nutrition Parameters: N Concentration versus N Accumulation

The estimation accuracies for LNC and PNC were low when using VIs and the best-performing VI could explain no more than 30% of the variability in N concentration. However, the estimation accuracies for LNA and PNA were significantly higher than those for LNC and PNC. The optimal VI (NDRE) could explain 77% and 73% of the variability in LNA and PNA, respectively. The LNC (or PNC) varied in a narrow range across the growing season of rice crops. It decreased from the beginning, leveled off in the middle and then decreased until harvesting [6]. This trend could be characterized by NDRE with a linear but weak relationship. However, LNA (or PNA) varied in a relatively wider range and kept increasing across the entire season. Particularly, NDRE tended to saturate when LNA (PNA) increased to a certain level. A linear regression would not be able to capture this decreasing rate of nitrogen accumulation. Therefore, LNC or PNC exhibited weak relationships with VIs in linear regressions and LNA or PNA in exponential regression [3,5]. Besides, the N absorption features are located in the shortwave infrared (SWIR) region rather than the visible and NIR region [30]. Li et al. [31] improved crop LNC estimation from the SWIR reflectance spectra of fresh leaves through enhancing the absorption features of nitrogen in the SWIR region. However, those VIs derived from UAV multispectral imagery were all based on visible and NIR bands, which could be easily affected by chlorophyll content and canopy structure [32]. Furthermore, the N dilution effect might be another reason for the low estimation accuracy of LNC and PNC with VIs due to the decrease of N concentration along with the increase of biomass [33]. The good relationships between NA and VIs might be attributed to the strong capabilities of VIs in retrieving biomass [34] and the close correlation between N accumulation and biomass (Table 3).

In texture analysis, the optimal NDTI for LNC and PNC estimation were all composed of NIR and RE bands with similar texture metrics (e.g., CON, DIS, COR). Although both LNC and PNC could be estimated with a common index NDTI (DIS_{B4_D3}, COR_{B5_D3}), the relationships between texture indices and NC were not strong in both linear and exponential regressions. That might be because leaf color turned light to dark repeatedly in response to N fertilization at the vegetative stages and kept yellowing at reproductive stages [6]. This reduced the sensitivity of texture metrics to the heterogeneity of tonal variation caused by N status in crops across the whole growing season. Specifically, NDTI(MEA_B5, MEA_B4) and NDTI(MEA_B5, MEA_B2) were the optimal texture indices across the whole growing season in exponential regression for LNA and PNA estimation, respectively (Table 5). Because allometric variation of LNA or PNA in the season was largely attributed to biomass (Table 3), texture indices could characterize the biomass related variation with an exponential regression [7,8].

Furthermore, the variables in the multivariate model for LNA included one N-sensitive variable (NDTI(DIS_{B4_D3}, COR_{B5_D3})) and two biomass-sensitive variables (NDRE and NGRDI). Therefore, LNA could reflect part of the variation in LNC [35,36], but the exact proportion needs to be investigated in future work. However, the optimal multivariate PNA estimation model consisted of entirely AGB-sensitive indices such as NDTI(MEA_{B5_D3}, MEA_{B2_D4}), MTVI2 [29], and NGRDI [20]. This also suggests that the variation in PNA might be dominated by biomass and contain weak information on N concentration which corresponds well with the findings in Reference [5].

4.2. Directional Effect of Texture Analysis on Row-Planted Crops

The directional effect of texture analysis was rarely investigated in the existing literature, since the majority of previous studies executed texture analysis with the default direction (45°) [7,15]. Some studies calculated texture metrics with different directions but did not explicitly explain the reason [10,16]. That might be because most of them studied naturally grown forests, and the trees were distributed disorderly. In contrast, we found that texture metrics had a significant directional effect (Figure 2). That is because rice plants are grown in rows, and the local window sliding along the row orientation contains more homogeneous vegetation than from other directions. Furthermore, the rice plants in the same row grow more homogenously than those in different rows, resulting in lower contrast and higher correlation in the along-row direction than in other directions (Figure 5). Texture indices calculated with texture metrics from different directions had different performances on the N nutrition parameter estimation (Table 6). However, the texture feature MEA represents the average values within the moving window and the textural information extracted from the UAV images is the average value from an ROI in the non-sampling area. As a result, MEA was not affected by directions, and the same R² value was obtained with the optimal NDTIs composed of MEA features calculated in different directions for LNA and PNA estimation (Table 6).

Texture metrics calculated in directions D3 and D4 constructed the optimal NDTIs for LNC and PNC estimation which might be explained by the strong capability of across-row texture metrics in differentiating the tonal variations caused by N status. Moreover, the estimation accuracy (RMSE = 0.24, RRMSE = 10.52%) of LNC decreased marginally if the texture indices in the multivariate model of LNC were derived from the direction D3 alone. Additionally, the same estimation accuracies were obtained for other N parameters (LNA, PNC, and PNA) through the multivariate models with all the texture indices calculated in the direction D3 when compared to the original multivariate models (data not shown). Therefore, all N nutrition parameters could be estimated at nearly the highest accuracies with the combination of VIs and texture indices calculated with texture metrics in the direction D3 which could simplify the use of texture analysis significantly.

Texture analysis also involved the optimal selection of texture calculation algorithm, spectral band, and window size [14]. Although window size has a considerable effect on the estimation accuracy of forest biomass due to the mismatch of spatial scale between remotely sensed pixel size and tree canopy [14,16], it showed no significant influence in the present study. This was because most window sizes contained a large proportion of crops due to the large canopy coverage since the jointing stage. Rice crops are often planted in the row distance of 24–30 cm, and multispectral images collected at 100 m usually possess a spatial resolution of 5–7 cm. Hence, a larger window size might be a good choice for texture analysis at early growth stages. Yue et al. [7] also found the optimal image resolution for using image textures to estimate AGB in winter wheat depended on the crop canopy size and row spacing. Therefore, the optimal window size has to be taken into consideration in terms of row space and image spatial resolution for other row-planted crops (e.g., soybean, corn).

4.3. The Benefits of Fused Information for Enhancing N Signals

Numerous studies were dedicated to improving N concentration estimation with different approaches due to the importance of N nutrition status in crop management. The most commonly used way was to propose new VIs for N concentration estimation [37,38]. Although those VIs yielded high estimation accuracies as shown in the literature, the unstable performance was reported in other studies [39]. In this study, VIs had poor performance in N concentration estimation, which is in line with References [3,5]. Texture indices showed comparable performance when compared to VIs which is in contrast to the performance of texture ratios in forest AGB estimation [10,16]. However, a significant improvement in LNC estimation was obtained when using the combination of texture indices and VIs as compared to using VIs or texture indices alone, with an increase of more than 35% in RMSE (Table 7). This finding agrees well with the results of References [7] and [8] which reported that the combination of texture indices and VIs improved the estimation of crop AGB significantly. In contrast to relevant studies [6,40], a universal model could be used to estimate N concentration across the entire season, which was a substantial improvement in crop N status monitoring with remotely sensed data. Compared with the new spectral index proposed by Stroppiana et al. [38] with hyperspectral data (R² = 0.65), our multivariate model derived from multispectral data could even yield a similar R² value (R² = 0.68) for PNC. A universal model suitable for the entire season could not only be used to guide N fertilization applications at the early growth stages [1,41] but also to predict crop yield and quality before harvest [42,43].

The improvement in LNC and PNC estimation induced by the addition of texture indices stemmed from the enhancement of N signals for the early stages. VIs could not be used to estimate N concentration at the early growth stages (from tillering to booting) with a universal model, because N signals were hampered by the rapid increase of biomass and the large proportion of background materials [5,6]. However, leaf color changed remarkably due to the N fertilizer applications conducted at jointing stage [11,12], and texture metrics could characterize the spatial distribution of tonal variations caused by N status. The values of CON feature had a greater variation between jointing and booting stages than reflectance (Figure 6). Furthermore, VIs had weak capability in N status monitoring at reproductive stages due to the high canopy coverage [2]. However, leaf tone had a visual color change from the vegetative period to the reproductive period, and the values of CON and DIS were significantly different at reproductive stages (Figure 6). NDTI has a wider variation than NDRE at high N levels at late growth stages (Figure 7). Therefore, the N signal could also be detected by those texture metrics at the late growth stages.

Moreover, texture metrics have the capability of smoothing the spatial heterogeneity between vegetation and background materials with a sliding window [9,16]. Therefore, the complementary information between texture and VIs might be useful for solving the saturation problem and reducing background interference. Although a similar conclusion has been drawn by Yue et al. [7] and Zheng et al. [8] for AGB estimation, this study represents the first solid evidence on improving the estimation of N nutrition parameters, especially leaf and plant N concentrations.

Multi-source data fusion was commonly used in previous studies to improve the estimation of crop biomass [44,45,46]. For example, Yue et al. [45] improved wheat biomass estimation by combining VIs and plant height derived from crop surface models (CSMs), but the construction of CSMs was tedious and time consuming. Although the fusion of data from different platforms could also acquire higher estimation accuracy [44], the data (e.g., LiDAR data) processing needs more professional skill. However, UAV-based image texture and spectral data are more convenient to acquire over small fields. Furthermore, the fusion of texture and spectral data could improve the accuracy of N status estimation of multiple growth stages, which solves the problem that N status was often estimated at only one growth stage or a short time window [6,40]. Therefore, monitoring N status across the whole season could be realized through the integration of texture and spectral information into the crop growth monitoring systems.

There might be some concerns that it is challenging to build a universal model for predicting the N status of rice plants across the entire growing season [3]. With spectral information alone, VIs often exhibit weak capabilities of detecting N status at late growth stages due to the fact of crop senescence. Textural information could capture the variation of leaf color within and between plots, especially at the late stages when leaf color changes dramatically. However, our results demonstrated that the fusion of spectral and textural information could improve N status monitoring efficiently. Furthermore, the multivariate models for LNC and PNC were tested on different dataset groups with satisfactory validation performance (Table 8). The applicability of those models still needs to be improved through further testing with more datasets from different geographic sites.

4.4. Potentials for Other Platforms

In this study, we found the combination of spectral and texture information in NIR and RE bands was superior to that in other bands in rice N status monitoring. As for monitoring crop N status in large areas, satellite imagery has been widely used but with low estimation accuracy. That is because satellite imagery has low spatial resolution and low spectral resolution normally with blue, green, red and NIR wavebands. RapidEye is the first launched satellite with a RE band and successfully applied to many aspects of precision agriculture [47,48]. Other satellites with RE bands (e.g., Sentinel-2, WorldView-2, and Gaofen-6) have been launched and have yielded significant improvement in the estimation accuracy on agronomic variables [49,50]. These findings from satellite observations were consistent with our results that the RE-based VIs were superior to other indices.

To date, texture analysis on satellite imagery was mostly used for forest AGB [14,15,16]. Our findings provide a strong reference on crop N status with texture analysis of satellite images equipped with the RE band. However, in terms of spatial resolution and crop row spacing, the spatial resolution of satellite imagery (highest spatial resolution is 0.3 m by WorldView 3) is insufficient for distinguishing the narrow row spacing in cereal crops (e.g., rice, wheat). Furthermore, crops could be planted in multiple directions over large areas and the direction effect of texture metrics may not be significant in satellite images with lower spatial resolutions. Given the suitability of UAV platforms for small-scale applications for the moment, the direction of texture analysis could still exist in UAV imagery. When using texture analysis on row-planted vegetation (e.g., staple crops, vegetable crops, and fruit trees) at the farm level with high resolution images, it is still beneficial to take consideration of the direction effect for improved nitrogen nutrition monitoring.

5. Conclusions

We investigated the potential of texture analysis from UAV-based multispectral imagery for N status monitoring in terms of texture metrics, texture directions, and texture indices. For N concentration estimation, low estimation accuracies were obtained using VIs, individual texture metrics or texture indices alone. However, considerable improvements were achieved with the combination of VIs and texture indices. The multivariate models with fused information yielded the highest accuracy for LNC (RMSE = 0.22 and RRMSE = 9.70%) and PNC (RMSE = 0.19 and RRMSE = 15.34%). Moreover, the multivariate models exhibited high estimation accuracy when tested on dataset in different years and growth stages. For N accumulation estimation, NDRE yielded high estimation accuracies on LNA (RMSE = 1.35 g m⁻² and RRMSE = 24.56%) and PNA (RMSE = 3.16 g m⁻² and RRMSE = 32.20%). Compared with the models based on VIs or texture indices alone, the multivariate model also yielded higher accuracies for LNA (RMSE = 1.16 g m⁻² and RRMSE = 21.04%) and PNA (RMSE = 2.38 g m⁻² and RRMSE = 24.17%) estimation. When calculating the texture metrics, the across-row direction (D3) should be used for the best performance over row-planted crops instead of the default diagonal direction (D2). These findings could serve as useful references for deriving appropriate texture indices from multispectral UAV imagery for N nutrition monitoring. This study provides new insights on enhancing the N signals from rice canopies by adding the textural information to conventional VIs. It has great potential in conducting rapid, accurate monitoring of N status monitoring over other crops and developing effective practices of improved crop management.